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ABSTRACT 


A numerical algorithm was developed for computing the 
nonlinear dynamic response of a ring-stiffened, nearly cir- 
cular cylindrical shell of finite length under transient, 
axisymmetric radial@loads of arbitrary axial distribution 
Nonlinear Donnell-type equations were solved using Fourier 
series expansions of the dependent variables in the circum- 
ferential coordinate, modified finite difference approxi- 
mations of the axial derivatives, and Newmark's beta-method, 
combined with Gauss elimination, for the time integration. 

The response of a simply supported shell under an 
exponentially decaying, uniform pressure was computed for 
peak pressures and total impulses between the static 
buckling limit and the dynamic buckling limit. Near the 
dynamic buckling limit, the exponential growth of the static 
buckling modes dominated; but as the peak pressure was 
reduced, the parametrically excited Mathieu modes became 
increasingly important. The significance of Ganping, 9eme 
initial imperfections, and nonlinear coupling was also 


investigated. 
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CHAPRER I 
INGROBUCT ION 
I. SCOPE OF THE INVESTIGATION 


This dissertation considers the general nonlinear 
dynamic response of a nearly circular, nearly cylindrical, 
ring-stiffened shell of finite length subjected to a tran- 
sient, axisymmetric radial load. 

The research»was part of a project for the Naval 
Weapons Laboratory, Dahlgren, Virginia, to analytically 
predict the dynamic response of the large NWL shock 
tube a} under operational loading conditions. The various 
sections of the shock tube can be idealized as thin, ring- 
stiffened, circular cylindrical shells having either clamped 
or simply supported ends. The operating loads can be de- 
scribed as internal, propagating, axisymmetric pulses having 
pressures both higher and lower than ambient atmospheric 
pressure. In the unloaded condition, the midsurface of the 
shell deviates from a perfect circular cylinder, but the 
unloaded shell was assumed to be stress free. Due to this 
initial imperfection of the shell, the dynamic response will 
not be axisymmetric even though the loading may be axi- 


symmetric. 


lNumbers within brackets refer to items in the List 
of References, page 179. 


2 


Although there has been a great deal of recent 
interest in problems of this nature ,~ there are, in fact, 
no previous analyses that are fully capable of solving the 
problem. 

Consequently, an algorithm was developed for the 
numerical solution to the Donnell-type, nonlinear differen- 
tial equations that describe the general dynamic response 
of the idealized shock tube sections to a transient, axi- 
symmetric, radial load of arbitrary time history and axial 
distribution. Some of the features of the method were 
demonstrated by computing the dynamic response of a simply 
supported shell under an exponentially decaying, uniform 
pressure. The peak pressures and the total impulses of the 
loads considered fall between the static buckling limit and 
the dynamic buckling limit. In addition, the Signiticame— 
of viscous damping, of the initial imperfections, and of 


the nonlinear coupling was investigated. 


rhe analysis is applicable to structures other than 
the NWL shock tube. For instance, a missile or submarine 
engulfed by a blast front presents a similar structural 
dynamics problem. 
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Pie LoLOREG ALA REVIEW 


Rings and Unstiffened Cylinders 

The nonlinear dynamic response of circular cylindrical 
shells has received considerable attention over the past 
twenty years, with the Russians doing much of the early 
work in the field. For example, the 1956 book by the 
Russian V. V. Bolotin [2] is a classic in the field of 
dynamic instability. In it the works of most of the early 
RusSian investigators are described in some detail. 

One of the earliest articles to appear in the American 
literature was a 1961 translation from the Russian of a 1958 
paper by V. L. Agamirov and A. S. Vol'mir [3]. In this 
paper the Bubnov-Galerkin method was used to obtain an 
approximate solution for the response of a closed cylindri- 
cal shell subjected to a hydrostatic loading that increased 
linearly with time. They assumed a function for the spatial 
variation of both the initial imperfection and the total 
deflection, allowing for variation in the axial and circum- 
ferential directions. Thus, the problem was reduced to the 
numerical solution of a single nonlinear ordinary differen- 
tial equation in the time dependent amplitude of the dis- 
placement function. 

me 1 96a, Jo CrevYaow [4] published one of the @mHrst 
American papers to appear in the field. He investigated 
the stability of an infinite cylinder under a uniform radial 


pressure applied as a rectangular shaped function of time. 


2 


The motion was assumed to consist of a uniform radial dis- 
placement plus one flexural mode having two waves around 
the circumference. The only coupling retained was the 
influence of the axisymmetric mode on the flexural mode. 
The flexural mode was shown to respond in an oscillatory 
or exponential manner depending on the magnitude of the 
pressure. 

The Agamirov and Vol'mir paper [3] was followed in 
1962 by another translation from the Russian of a 1960 
paper by Y. I. Kadashevich and A. K. Pertsev [5]. They 
investigated the response of a finite cylindrical shell 
loaded by uniform radial pressure applied as a step func- 
tion in time. They used combinations of two axisymmetric 
modes and one nonaxisymmetric mode to represent the initial 
imperfection and the radial displacement. However, these 
modes did not satisfy the boundary conditions. Using 
Lagrange's equations, ordinary, coupled, nonlinear differ- 
ential equations were obtained which accounted for the 
radial inertia of all three oscillatory modes. These 
equations were integrated numerically to obtain the time 
histories of the modes. The critical load was defined as 
the load which caused stresses in excess of the ultimate 
strength of the shell material. This paper is considered 
most significant because it retained all the nonlinear 
coupling between the three modes. Most later works do not 


retain all of the coupling terms. 
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Another very Significant paper was published by 
Yao [6] in 1963. He investigated the dynamic stability 
of closed cylindrical shells under several combinations 
of harmonic and static loads. He assumed the radial 
motion to consist of a uniform quasi-static displacement 
plus a small, nonaxisymmetric, flexural perturbation that 


3 NS ities Carl ter 


Satisfied the boundary conditions. 
work [4], Yao retained only the influence of the axi- 
symmetric mode on the flexural mode. The governing equa- 
tion for the motion of the perturbation was an ordinary 
differential equation with a harmonically varying coeffi- 
cient, namely, the classical Mathieu equation.” When a 
flexural mode corresponded to a point in an unstable region 
of the Mathieu stability chart, the response of the shell 
was assumed to be unstable. Yao concluded that a shell not 
only can withstand static loads up to the classical 
buckling load, but also can endure an additional harmonic 


load so long as the harmonic component does not render any 


of the flexural modes unstable. 


3The radial inertia of the axisymmetric mode was 
neglected. 


‘the solution to the Mathieu equation is either 
oscillatory and bounded, or oscillatory and unbounded, 
depending on the parameters [7]. The unbounded condition 
1s known aS parametric resonance or dynamic instability. 
The distinction between a bounded response and one that 
is unbounded is determined from the Mathieu stability 
Chane. 
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In late 1963, J. D. Wood and Lhe wigs. ae lem 
sented the results of their investigation of the stability 
of an axially loaded shell subjected to a time dependent 
uniform radial pressure. Their assumptions and method 
of solution were essentially the same as those employed 
by Yao [4, 6]. However, they did not neglect the radial 
inertia of the axisymmetric mode as Yao had done. 

In 1964, J. N. Goodier and I. K. McIvor [9] published 
a paper which considered the response of an infinitely 
long shell subjected to nearly axisymmetric, radial impulse 
of low intensity. Mathieu's equation governed the early 
growth of the flexural modes. Since the impulse was of 
low intensity, a Single dominant flexural mode could be 
identified from the Mathieu stability chart. They assumed 
the radial displacement to consist of a uniform radial 
mode plus the single dominant flexural mode. The fully 
coupled motion of these two modes was computed numerically. 
During the motion, the axisymmetric mode and the flexural 
mode exchanged energy in a cyclic manner because of the 
total coupling between the two modes. Although the 
flexural mode corresponded to an unstable point on the 
stability chart, its motion remained bounded since only 
a finite amount of aHewey was imparted to the shell by 
the load. This paper demonstrated that the motion of a 
flexural niece 1s not necessarily unbounded even though 
its early growth is associated with a point in an unstable 


region of the Mathieu stability chart. 
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H. E. Bindberge [10], im a companion to the abewe 
paper, investigated the dynamic buckling of thin shells 
under very intensive impulsive loads. He found that many 
short wavelength flexural modes were excited to large 
amplitudes, and that buckling occurred as the result of 
hyperbolic growth, rather than any oscillatory growth. 

In 1964, R. S. Roth and J. M. Klosner [ll] investi- 
gated the dynamic response of a cylindrical shell subjected 
to an axisymmetric, axial step load. They used a Donnell- 
type shell theory and assumed that the radial displacements 
and the initial imperfection were comprised of a linear 
combination of four discrete modes, all of which did not 
satisfy the boundary conditions. Radial inertial forces 
in all modes were considered. The Galerkin-Ritz method 
was used to obtain the coupled, ordinary, nonlinear 
differential equations in the time dependent coefficients. 
The resulting equations were then integrated numerically. 
Prior to this work, the Russians Kadashevich and Pertsev [5] 
were apparently the only ones who had accounted for the 
Simultaneous, coupled motion of more than one flexural mode. 

In 1965, a survey article by R. M. Evan-Iwanowski [12] 
reviewed the state of the art on the parametric resonance 
of structures. The bibliography is very comprehensive with 
numerous references to the early Russian work. In this 
article Evan-Iwanowski suggested that future work should be 
devoted to study of the interaction between the dynamic and 


parametric responses. 
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M. P. Bienick, T. C. Fan, and L. M. Lackman's 1966 
paper [13] investigatedmthe dynamiewstabilieot & sei 
drical shell under a uniform radial pressure, applied 
harmonically or as a ramp function of time. The axisym- 
metric deflection state varied as a half sine wave along 
the axis, thus satisfying the simply supported boundary 
conditions. One flexural mode of arbitrarily small magni- 
tude was assumed to comprise the perturbed state. It 
also satisfied the boundary conditions. The effects of 
the radial inertia in both modes were included. The 
stability problem was reduced to Mathieu's equation for 
both types of loads. Under the ramp loading, they found 
that not one, but several flexural modes were unstable 
if the shell was short. 

The first works devoted to the dynamic stability of 
a cylindrical shell loaded by a moving axisymmetric load 
were by G. A. Hegemeir [14, 15] in 1966 and 1967. The 
Shell was infinite in length, and the load moved at con- 
stant velocity, leading to a steady state condition. 
Hegemeir, like Yao [6] and Bienick, Fan, and Lackman [13], 
investigated the stability of the nonaxisymmetric response 
with respect to arbitrarily small perturbations. However, 
the governing equation was not the Mathieu equation. 

In 1966, T. C. Fan [16] retained the axial and circum- 
Par cenin inertia in an investigation of the dynamic 
Stability of a simply supported cylinder loaded by uniform 


radial pressure applied as a ramp function of time. His 
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results showed that neglecting the axial and circumferential 
inertial forces has little effect on the stability of small 
flexural perturbations. This seems to be the only study 

of the significance of the axial and circumferential iner- 
tial forces. 

In 1966, M. H. Lock [17], like Goodier and McIvor [9], 
studied the response of a ring, or an infinitely long 
cylinder, but under step loads instead of impulsive loads. 
Using finite deflection theory, he computed the coupled 
motion of the axisymmetric mode and one nonaxisymmetric 
mode. The conditions for instability to small perturba- 
tions were shown to be necessary for the existence of large 
displacements, but not sufficient. Significant to the 
present investigation, he showed that a flexural mode could 
be excited to large amplitude in either of two ways: (1) 
When the magnitude of the load exceeded the static buckling 
lead of “a mode, that mode grew in essentially an expo- 
nential manner with time until limited by nonlinear effects. 
(2) A flexural mode could also be excited to large ampli- 
tudes by the oscillation of the axisymmetric mode, and the 
growth was an oscillation of increasing amplitude, as in 
parametric resonance. 

A Pweowrepert by H.wEy hindberg, ee al. [ie and ‘em 
sulpequenten 1968) article by D. L. Anderson and H. E. 

Lindberg [19] are of considerable interest to this disser- 
tation. They presented two theoretical models, a tangent 


modulus model for plastic behavior and an elastic model, 
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for predicting the dynamic loads required to buckle a 
cylindrical shell. Only the elastic model is of interest 
here. This model made use of a Donnell-type shell theory 
which accounted for initial imperfections. The radial 
motion was assumed to consist of a uniform radial displace- 
ment plus nonaxisymmetric flexural modes that satisfied 

the simply supported boundary conditions. The spatial 
distribution of the initial imperfection was taken to be 
the same as that assumed for the flexural mode shapes. 
Following Yao [4, 6], Wood and Koval [8], and Bienick, 

Fan, and Lackman [13], the influence of the axisymmetric 
mode on the flexural modes was the only coupling retained. 
No coupling was retained which might extract energy from 
the uniform radial mode, nor was any coupling between 
flexural modes retained. An exponentially decaying uniform 
load and a linearly decaying uniform load were con- 
Sidered. The maximum displacement of several flexural 
modes was found by numerical integration, and buckling 

was assumed to have occurred whenever the maximum displace- 
ment of any mode exceeded one thousand times the initial 
imperfection in that mode. Only the lower frequency, 
exponentially growing buckling modes were considered. The 
parametrically excited modes were not treated. The authors 
reasoned that these modes would dissipate energy locally by 
plastic flow and, hence, would eventually dampen out before 


causing deflections large enough to buckle the shell. 
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The most recent work is the paper by I. K. McIvor and 
E. G. Lovell [20] which appeared in late 1968. They 
studied the dynamic response of finite length cylindrical 
shells subjected to a uniform radial impulse. They assumed 
a set of modal displacement functions that simplified some 
of the algebra, but did not satisfy end conditions of 
usual interest. By first using a small nonaxisymmetric 
perturbation approach that retained only the influence of 
the axisymmetric mode on the flexural mode, Mathieu's 
equation was obtained. They then determined from the 
Mathieu stability diagram the four or five flexural modes 
having the highest initial growth rate. Using these 
modes and the axisymmetric mode, and retaining all 
coupling between modes, they employed the Runge-Kutta 
integration method to compute the response of each mode 
to the uniform radial impulse. A slight irregularity 
in the impulse was used to initiate the response of each 
nonaxisymmetric mode. They found that the modes under- 
went beat-like oscillations accompanied by a continual 
exchange of energy between modes. The stress levels 
attained during the nonsymmetric motion were much higher 
than those developed in an axisymmetric dynamic response 


to the same load. 


Dynamic Stability Ok Ringacite nkened Shells 


The treatment of discrete ring stiffeners in a non- 


linear dynamic analysis would provide a suitable topic for 


3/3 


a dissertation in itself. In the course of this research, 
a brief study of the problem was conducted, and a seemingly 
Suitable method of analysis was selected. It is not the 
only method, and perhaps not even the best method. A 
complete review of the related literature is deemed to be 
beyond the scope of this study. Only a few references 

are mentioned to provide a starting point for an interested 
reader. 

In a 1967 dissertation, W. K. Dietz [21] investigated 
the dynamic response of eccentrically reinforced, circular 
cylindrical shells subjected to dynamic axial loads. He 
considered only shells with closely spaced stiffeners 
and treated the structure as an equivalent orthotropic 
shell. The Galerkin method was used to obtain a set of 
ordinary, nonlinear differential equations in time which 
were integrated by employing a predictor-corrector method. 

J. M. Klosner and H. N. Franklin [22] presented a 
report in 1968 which dealt briefly with the dynamic insta- 
bility of long reinforced cylindrical shells under axial 
loads. They, like Dietz, treated an equivalent ortho- 
tropic shell and used the Galerkin method to obtain the 
governing differential equations of motion. The equations 
were integrated by the Runge-Kutta method. 

The analysis of the rings contained in this disser- 
tation is essentially an adaptation of a method used by 


D. L. Block [23] in a 1968 report dealing with the influence 
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of discrete ring stiffeners and prebuckling deformations 


on the static buckling of axially loaded shells. 
TIT. NEW FEATURES OF THE INVESTIGATION 


The studies which have been reviewed in the proceding 
section appear inadequate for one or more of the following 
reasons: 

1. The motion can not be predicted after onset 
of nonaxisymmetric motion. 

2. The boundary conditions are not satisfied 
by the assumed displacement functions. 

3. Coupling due to nonlinear effects is only 
partially retained. 

4. The assumed displacements are not general 
enough; or an insufficient number of modes 
have been included to encompass both the 
higher frequency, parametrically excited 
modes and the exponentially growing, lower 
frequency, buckling modes. 

5. The analyses do not handle loads which depend 
on the axial coordinate. 

The algorithm developed in this dissertation is more 
general than any of the earlier analyses and does not con- 
tain any of the above shortcomings. This is accomplished 


by the use of the following new features: 


oS) 


ie 
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The displacements are expanded in Fourier 


series in the circumferential coordinate. 
The series are then truncated, and only 
the significant modes are retained. 

Those modes which do not respond can be 
discarded. As many as fifteen modes are 
included, and this number can be increased 
at the expense of increased computing time. 
In this study all of the nonlinear coupling 
between the modes is retained, and enough 
modes are included to simultaneously en- 
compass the parametrically excited modes 
and the exponentially growing buckling 
modes. 

set of modified finite differences is 
employed for the derivates with respect to 
the axial coordinate for the first time 

in a nonlinear shell problem. These 
modified differences are shown to be 
superior to the conventional difference 


approximations to derivatives. 


The iterative Newmark beta-method is employed 


for the time integration of the radial 
equations of motion. The axial and circum- 
ferential displacements are obtained by 


Gauss elimination during each iteration. 
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4. The boundary conditions are consistently 
observed. 

5. Axisymmetric radial loads of arbitrary time 
variation and axial form can be readily 
handled. 

6. Initial imperfections symmetric about a 
plane through the axis and of arbitrary 
axial distribution are accounted for. 

7. Viscous damping forces are included in the 
equations of motion. 

8. The effects of discrete, eccentrically 
located ring stiffeners have been included 


in the analysis. 


IV. ORGANIZATION AND PREVIEW 


The main body of the dissertation is divided into six 
chapters. The present introductory chapter is followed by 
the development of the equations of motion in Chapter II. 

In Chapter III, a numerical model is developed that replaces 
the partial differential equations of Chapter II with 
difference-differential equations in the time variable. In 
Chapter IV, methods for generating a solution to the 
difference-differential equations are considered, the com- 
puter program is briefly described, and the tests which 

were conducted to verify the computer program are discussed. 
Numerical results for a problem, a problem similar to one 


considered earlier by Anderson and Lindberg [19] and McIvor 


a7 


and Lovell [20], are presented in Chapter V along with 
studies of the significance of damping, of the initial 
imperfection, and of nonlinear coupling. Chapter VI 
contains a summary of the work, the new findings, and 
recommendations for future work in the field. 

Four appendices contain some of the minor details 
of the development of the algorithm and a complete docu- 
mentation of the computer program, including instructions 
for preparation of input data cards and other information 


necessary for the operation of the program. 
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CHAPTER II 


PORMULATION OF THE PROBLEM 


I. INTRODUCTION 


In this chapter the governing nonlinear differential 
equations for the dynamic response of a ring-stiffened, 
nearly circular cylindrical shell of finite length are 
derived. A Donnell-type nonlinear theory [24, 25] that 
accounts for initial imperfections in the shape of the shell 
is used to describe the shell response; the discrete 
circular rings are treated by linear beam and elementary 
torsion theories. The shell and ring materials are 
assumed to be homogeneous, isotropic, and elastic; and the 
rings may be fabricated from a different material than the 
shell. The centroids of the ring cross sections are eccen- 
trically located with respect to the midsurface of the 
Shell. Radial inertial forces of the shell and the rings 
have been included; but their axial, circumferential, and 
rotatory inertial forces have been neglected. Viscous 
type damping has been included to account for dissipative 


Cnet en 


its COORDINATE SYSTEM 


The motions of the shell of length L will be referred 
to the fixed coordinate system shown in Figure 2.1. The 


positive direction of the axial coordinate x, the 
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Midsurface of 
Perfect Shell 


FIGURE 2. 


GEOMETRY AND COORDINATES 
OF SHELL 
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circumferential coordinates y and 9, the radial coordinate 
z, the axial displacement U, the circumferential displace- 
ment V, and the radial displacement W* are as indicated in 
Figure 2.1. The origin of the coordinate system is located 
at the midsurface of a perfectly circular cylinder of 
radius a. The ends of the shell are located at x = 0 and 

x = L. The shell thickness is denoted by h. 

The motion of a ring will be described by the dis- 
placements of the centroid of the ring cross section ua 
yt and Wee and the rotation about a normal to the cross 
section 8, with positive directions as indicated in Figure 
2.2. The distance e between the centroid of a ring cross 
section and the midsurface of the shell is positive when 


the rings are on the outside of the shell. 
IIIT. STRAIN-DISPLACEMENT RELATIONS 


Shell 
In the Donnell-type nonlinear theory, the relations 


between the midsurface strains and displacements are given 


oh easy: 
i Z 
ene to) (WN |) er Wy We, (Zena) 
x x xX 
W 2 — 
Se, eee ee CW, + W, WwW, (22D) 
y y = ear e y ¥ 
a = UTS + Mins + Wie Wry + We Wiy + NO (22e) 
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FIGURE 2.2 
DISPLACEMENTS AND ROTATION 
OF A TYPICAL RING 
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X, = W, G2 eed.) 
ae = Ny (2.le) 
X = W, CZeP te ) 


where €, and Ee, are the extensional strains, a oe is the 
shear strain, or and x are the curvature changes, Xy¥ is the 
twist, W is the radial displacement of the midsurface of 
the nearly circular cylinder, and W is the radial deviation 
of the unstrained midsurface from the midsurface of a 
perfectly circular cylindrical shell. Hence 

we =Wtw 


Here, and in the sequel 





; = yx 
ae 
‘yy 2 
oy 


Rings 

The rings are assumed to be circular in shape and 
radially thin compared to their radius a”. The strain- 
displacement relations taken for the rings are the ones 
used by G. D. Galletly [26] in his investigation of the 
free linear vibration of ring-stiffened shells. The 


strains of the centroid of the ring are given by 
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1 il ’ 

X1 = aa (8 + oe. 
a a 

Xp = ny (WE + Whog) 
(a) 

(223 
r 

Ba lee 
at “6 at 

ae ale ic ceecr eis 

Eo = at (Vig W ) 


where X7 is the curvature change in the plane normal to the 
ie TeGly X5 is the curvature change in the plane of the ring, 
@ is the angle of twist per unit circumferential length, 
and E5 is the circumferential strain. 

Assuming 8 to be small, e<<a, and that W and Wry ange 


zero at the point of attachment, the displacements and 


rotation of the ring can be geometrically related to the 


shell displacements [26]. This leads to 
w= W 
B= Wee 
(229 
uv =U+e B 


vi=avVieu, 


and the radius of the ring 1s given by 
a” = ate 
Thus, the ring strains can be expressed in terms of the 
shell displacements. Substituting (2.3) into (2.2) and 
g 


NOLIng, Ehiate Je. — a ( lie gives 
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r ai a = 
= —-— oe + 
a ma eee “7 (Ur + © Wryyy)] 
r aE 2 
pee W, 
Xo 2 ( yy? 
(224) 
s Uae & Wry 
= ee + 
p “x ( Wry aaa ) 
a a ol 
EG “x (Vig + e Wivy =) 


Note that (224) does not contain Wy, the initial imper- 
fection of the shell. This is due to the fact that the 
intial imperfection W and its slope We, have been taken 


equal to zero where the rings are attached. 
Ey. VCCNSTITUTIVE RELATIONS 


The resultant forces and moments per unit length of 
shell midsurface are obtained by integrating the stresses 
over the thickness. The normal forces are N.. and oak the 
shear forces are Ny and ee the bending moments are M. 


and a and the twisting moments are My and M The 


ve 
positive directions of these forces and moments and the 
applied pressure P are as shown in Figure 2.3. 


The stress resultants are related to the strains by 


the relationships 


No = B(e, + er, (2.5a) 

Ny = B(Ey + View) (2.56) 
aa oe 

N = N =e 2 aoe 

ees een, ( 
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FIGURE 2.3 


POSITIVE DIRECTIONS OF FORCES AND MOMENTS 
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M = -D(X + VX) (2. 5a) 


x Xx 

=> + \V: 2.5e 
A = -D(x, Xx) ( 
Mey = Mx = (a) DX Sy (2.52) 


Maere V is Poisson's ratio, and B and D are the extensional 


and flexural rigidities, respectively, defined by 





ie: En, 
i. 
(2.6) 
, Eh? 
D = ——--- 
12 (l-v"} J 


Young's modulus is denoted by E. 
VY. E@UATIONS OF MOTION 


The governing equations of motion are derived in this 
section by means of Hamilton's Principle. For the case 
where nonconservative forces are present, Hamilton's 


principle states that [27] 
6 6 7 (OP+97). = 0 (2.7) 


where 2Pis the total work function, F is the total kinetic 


energy, and t is the time coordinate, ty and t.. being arbi- 


2 
trary times. 

The total work function is the sum of the work function 
of the external pressure and damping force Le the work 


fometion of the shell at and the work function of the 


rangs “oie. 
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Thus 
B= Yo + PP + #2 es) 


The variation of Lie US Galv enwby 


L 2Ta 


69° =-fsf sf (CW,, - P) 6W dy ax (2.9) 
O O 


ic 


where C is the viscous damping coefficient. The termscw, 


G 
and P are held constant during the variation. 
The variation of Ae 1S given by 
9° L 27a 
6 = - N ide. +N _ 6 + N 6¢€ 
se Oy ROR * NS Cy ee 
(2.10) 


= 1 == 
M,OXy oe ey ol oie dydx 


The in-plane forces and moments are also held constant 
during the variation. 

The contribution ae of the ring forces to the total 
work function is taken as the negative of the strain energy 


in the rings. Hence, 6 | jswgivenwby [26 


r 


1 Nr 21a r 2 r 2 
SW = - we EOS {ETL (4) + EPR 
j=1 O 
(2.58) 
2 r 2 = 
+ G Cop Cd )i ack EAL. (€,) ] dy p eee ax 
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where 


characteristic length 

moments of inertia of the ring cross section 
with respect to axes passing through the 
centroid of the ring 

cross sectional area of the ring 

St. Venant's torsion constant for the ring [28] 
Dirac delta function, not to be confused with 
the symbol used to denote a variation 

number of rings 

Young's modulus of the ring material 

shear modulus of the ring material 

mass density of the ring material 

subscript which associates a variable with 


the 4th ring; for example, Be is the 


x-coordinate of the 5th sing 
= ie 
dy = 2_ dy 
a 


The Dirac delta function has been introduced here to allow 


the integration of 5m. to be carried out over the length 


of the shell. 


This technique has been used by Block [23] 


and is necessary in order to combine MW, RR, and Be. in 


a double integration form. It requires the introduction 


of a characteristic length uw that corresponds to the width 


of a ring. 


The total kinetic energy of the structure is given by 


OF _ oe” 
V¥-Ys* vr — 
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where 7. is the kinetic energy of the shell and 7, is the 
kinetic energy of the ring. Only the radial velocities will 
be included in the kinetic energy of the shell and of the 
ring; the axial, circumferential, and rotatory inertia of 
both the rings and the shell are assumed to be negligible. 
Thus, the kinetic energy of the shell is given by 


oh Lay 2a 2 
G_ =, i (W,,)° dydx (2.13) 


where p is the mass density of the shell material. The 


kinetic energy of the rings is given by 





ee d 2 
ry) ee (2 5 


Substituting (2..eem¥and (2°92) “inte (2.7) yields 


to 
ay Gan 
Saf PM nity at Hg Le + se, amaete 


il 
oO 


(2a) 


Substymtuting (2.1). dmiteOue( Zee 0 )icmeale Atel Tats O gale) leliaeeeiia al 
(2 2 Digi? 10) » a2 dil, af2cli2) wanda. 14)_1n bom aS) aoe 
performing the variation yield the following equations of 


dynamic equilibrium 





eee 
N + N > W, au, 
1. Siey rye ey ly a peetllennry yyyy 
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VI. BOUNDARY CONDITIONS 


The variational procedure also yields the admissible 


boundary conditions [29]. 


Simply Supported 


The boundary conditions that correspond to simply 


Supported ends are 


DL 


w = 0 
(2. 095) 
ne 0 
N = 0 at x = 0,L 
x 
Clamped 
The boundary conditions that correspond to clamped 
ends are 
U = 0 
V = 0 
(2320) 
w = 0 
W, = 0 at rer L 
og 


VII. NONDIMENSIONAL EQUATIONS 


Nondimensional Variables 


The nondimensional variables and parameters that will 


be used in the analysis are defined by 


awe - i =e 
f= 3 : a 8 a 
ie 4g 
ars ( 2 5 ) 
eae Ve") 
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ae We w =a5 W emi 
2 
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Nondimensional Equations of Motion 
Expressing the equations of motion (2.16), (2.17), 
and (2.18) in terms of the above non-dimensional quantities 


leads to the three equations 


(W, a 


ane ecucn wen A ro0 * Y'99e6 ) 
? j=l 


rege. 


= Urgg) a] oa = 0 
(2.22a) 


C3 (Wr 4 te ren aig = Vigg) Otere.) = 0 


(Ze225) 
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Ce ee a ae ee eee + Neg Tore. 
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The in-plane forces N , N_, and i‘ 
x 


y 


can be expressed as func- 


tions of the displacements U, V, and W by combining the 


constitutive relations, 


displacement relations 


(2.5a)sul2. 


(2.1 ay=nGZaace 


5c), and the strain- 


In nondimensional 


form 
= sl 2 re 1 W 
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Nondimensional Boundary EConarcions 


Simply Supported. The boundary conditions for the 
simply supported shell are expressed by (2.19). The 
corresponding conditions on the nondimensional displace- 


ments v and w are 


v =e (2.24) 


w = 0 at €  —.070) (2525) 


The bending moment My is a function of the curvatures and 


was expressed by (2.5e) as 


My = = ID Tr VX (2.5e) 


ExpresSing the curvatures in terms of the displacements by 


using (2.1ld) and (2.le) gives 


M.. = BOY ne te eed 


At the ends of the shell W = 0 which implies that Nee = 0; 


hence, the condition that My O at x = 0 and x = L implies 


that Mis a = 0 at x = 0 and x = L. In nondimensional form 


= = 1 
aS 0 at € 0, (2.26) 


Since Ny = 0 at x = 0 and x = L, it follows that ny = 0 at 


—€ = 0 and € = 1. Hence, from (2.23a) 


ee 45 (Wg) + w,-W,, + VIV, o-wHs (wy 9) * + Wr Wie] = 0 


oe 


at €&= 0,1 (2 220) 
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Clamped. The clamped™boundery cond terense te. 0) mean 


also be expressed nondimensionally as 


u- == 0 (2.220) 
v = 0 (22.290 
w = 0 (2.30) 
Wee = 0 at £€ = 1051 (2335) 


Response Symmetric to Midspan. 9 When seinem load mace 


initial imperfections, and the response are symmetric with 
respect to the midspan of the shell, only the portion of the 
shell between € = 0 and & = 1/2 needs to be considered. The 
nondimensional conditions that express symmetry at the 


midspan are 


u = 0 (2.223 
Vig = 0 (2.388 
Wee = 0 (2.34) 
eee = 0 at §& = 1/72 (2353 


The three equations of motion, (2.22), together with 


the expressions for Nes Ngs and n (2.23), and the boun- 


EO. 
dary conditions constitute a formulation of the problem 
under investigation. The next two chapters are devoted to 


their solution. 
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CUAL R wot 
DEVELOPMENT OF THE NUMERICAL MODEL 
ree Uv rRODgel TON 


In this chapter the nondimensional, nonlinear govern- 
ing partial differential equations derived in Chapter II 
are replaced by a set of difference-differential equations 
in which time is the only remaining independent variable. 
The difference-differential equations are henceforth 
referred to as the numerical model. 

The model is developed by expanding the dependent 
variables in Fourier series in the circumferential coordi- 
nate. The Fourier coefficients in the series are functions 
of t and the axial coordinate €. Products of Fourier series 
arise wherever there are nonlinear terms in the governing 
equations. The expansion of these products to single series 
generates products of the Fourier coefficients that couple 
the modes. The derivatives of the Fourier coefficients 
with respect to the axial coordinate are then approximated 
by modified finite differences, leading to the set of 


difference-differential equations. 
It. FOURZER SERIES EXPANSIONS 


During the early phases of the present study, con- 
Sideration was given to a numerical model that approximated 


both the circumferential derivatives and the axial 
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derivatives by means of finite differences. To achieve good 


ee ea 


resolution of the strains and curvatures of the cylinder, 


the grid points must be close enough together so that there 
are several grid points covering the shortest wave length 
of the response. Based on the works of previous investi- 
gators [9, 10, 20], the parametrically excited modes are | 
expected to have many waves around the circumference. 
Consequently, a very large number of points is needed for 
adequate resolution of the parametrically excited modes. 
For example, a flexural mode with sixteen waves around 
the circumference is not unreasonable to expect. Taking 
the number of axial stations to be twenty and assuming that 
eight node points per wave length are needed around the ; 
circumference, the total number of node points is approxi- 
mately 2500. Satisfying axial, circumferential, and radial 
equilibrium at each node yields over 7500 simultaneous 
equations. The number of equations can be reduced by one 
half if the response is assumed to be symmetric with respect 
to some axial plane. Even so, the number of equations is 
still very large. 

However, if the displacements are expanded in Fourier 
series in the circumferential coordinate, instead of using 
finite differences, the circumEcrential Vabiaelvoneo ter 
normal modes of the cylinder can be represented exactly. 
Furthermore, those modes which do not respond can be ex- 
cluded. Thus, the number of equations can be kept small 


while retaining the. significant features of the coupling 


a6 


between modes. Expansions of the type employed here have 
been used previously by R. E. Ball [30] in a nonlinear 


analysis of statically loaded shells of revolution. 


Form of the Series Expansions 
The assumption is made that the initial imperfection 
of the cylinder can be expressed by the Fourier cosine 


Series 


w (&) cos né (eons) 


0 


w(&,@) = 
n 


lm 8 


This form somewhat limits the applicability of the analysis 
Since the initial imperfection of most shells cannot be 
expressed by the cosine series alone. However, the form 
is still sufficiently general to permit the retention of 
coupling effects that hitherto have not been retained. 
There are no conceptual problems involved in using the 
complete Fourier series. Only the cosine series is used 
here because of limitations imposed by the practical con- 
Siderations of computer storage space and computing time. 
In order to correspond with the form of w given by 
(3.1), the remaining dependent variables are expressed in 


trigonometric series as follows: 


Tee Gee) a), tes (fo) cose (3.2a) 
n=0 

Vie oer) = r ee (eee) Sn on 0 (327 bp 
n=L1 

w(&,8,T) = ae (—,T) cos née (3.2c) 
n=0 


a9 


n (£,0,t) = £ n” (&,t) cos né (32369) 
e n=0 E 
ae. | 
Ng (€,0,T) = 2 Ng CE Gimmie Smiio (3, 309 
n=0 
Mpg (619, T) = ain neg (E07) sin né (3.3c) 


When these expressions for the displacements and in- 
plane forces are substituted into (2.22) and (2.23), the 
linear terms can be grouped together in one series. However, ; 
the nonlinear terms give rise to products of series and 
cannot be grouped with the linear terms until they are 
expanded into single series. The products of series are 
of three types; products of (1) two cosine series, (2) a 
Sine series and a cosine series, and (3) two sine series. 
The details of the expansion of these types of products 
are described in Appendix B of Reference [30]. Employing 
(3.1), (3.2), and (3.3), the nonlinear terms contained in 


(2.22c) and (2.23) can be expressed by 


2 co 
ae = ee We cos 180) ( ‘ 2008 30)= aoe cos ng 
I= 
(3.4a) 
5 (W, i wee) eee 10) % : -swisin 4e)=2 cua cos n9@ 
@ i=l 4=1 n=0 
(3. 4i55 
WrrWrg = ( : wi Gos 72.6) { ; SS an JO8)= 2 Gees sin no 
izo §& j=1 n=1 
(3.4648 
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a co ae | co ; | 
W, i= We COS me al fT w,aeecos ‘eas ia r co cos n@ (3.4d) 
00 : 00 , 00 * 
Wr Wr a= | 5 -iw- sin i@)( 2 -jw Sin j@)= 2 Comp COS n@ 
i=l j=1 n=0 
(3.4e) 
Wrew, =( 2 Wr COS.01:0) (2 -iwe Sin j@)= = oe Sin né 
1=0 j=1 n=1 
(opt) 
Wr Wy = ( 5 ~iwt 51-71 sO al 5 we cos j8@)= 5 Crt Sin n@ 
6 i=l s=0 §& n=1 
(349g) 
( 


Serer ee) en, cos 0) 2 Wirggtwsge)cos 38] 


(324) 
an E cos née 
n=0 
n  (w,-,,tw,-,)=( I ae sin i6) [ 5 See) sin 36] 
een, SGun 3a Ee a Eee 
1 
(34a) 
aac 
=? nN cos né 
n=0 oe 
WN eee ee) = | 5 ae cos 16) [ 5 ~4* (wI4H) cos j@] 
Cae eee age = 
= = 
(3.49) 
= a cos né 
n=0 


The products containing the initial imperfection, or deri- 


vates of the initial imperfection, are not nonlinear since 
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w(&,8) is assumed to be known. However, the products con- 
taining w expand in the same way as the nonlinear products; 


hence, they are included in the above list. 
Reduction to Two Independent Variables 


Equations of Motion. Expressions for the Fourier 
coefficients of the in-plane forces ne, Ng: and Neg? can be 
obtained by substituting (3.2), (3.3) and (3.4a)-(3.4g) into 
(2.33), and equating coefficients of like values of n. 


Performing the operations yields 


To lenenan n Tdleen ies a2 i) n 
Be a ita seeping: w ) + Cyy + Chyx 
(3.5a) 
V(Cam + Chimp) eB, Lue 
n n n n n n 
= ~ + 
no nv Woe aa Com + Cire 
(3°. 5) 
n n a 
+ V(Cyy + Crxx) Nai Lee: , 
Lee ae n _ n n 
"E@ 2 Veg = eee Cry) 
(3, 563 


The equations of motion in the axial and circum- 
ferential directions can be expressed in terms of the 
Fourier coefficients by substituting (3.2) and (3.3) into 
(2.22a) and (2.22b), and by making use of (3.5). The 


Bes lows 
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lip ee nv,. Soe n¢un aa lla: Le n?ud + 
2 eer? #2} 4 
n n 
eae (en-1) Wee = errs (We - u a ea = 0 
i 00 ee eee (3.6) 
Can 1l+v n l-v n n 
z == —v,-, + A2 
De av 5 nUr¢ + 5 V EE d 
Ny 
2 n 2 on 
+ C,[n(en =) Wer n wy OteZe) = @ 
j=l i 
a ee ee (327) 
where 
ae =e tee folice . Guor 
OG, e Oe V (Com & Crer, 6? 
(30) 
oo 
he n n n & n 
2 n(Cyn + Ciym + Copy) 7 Wee 
and 
1 n n n 
2 n (Ch, aC ee DG Cee) | 
(359) 
l-v Nn ii n n 
oa + 
mae ‘aes Greco ieee xemeteaiteOY 


All of the nonlinear terms are contained in 1 and 2”. 
Examination of these two quantities reveals that the axial 
and circumferential equations are nonlinear in we only; 
there are no nonlinear terms which contain u’ or v”™. 
Byeasubsertutang (392c), (3.3b), and (3.4h)- (3.4 9)" inte 


(2.22c), the equation of radial equilibrium is obtained in 


the form 
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ce n 2 4 n nh | 
— - + = - + 
To * rere eo ee eee : 
9 
a C. Wes (Eb o e —aee 
j=1 2 J J & 0 
gi NT) n ob nN ‘ 
+n.+n +n. + Ez R. 6(€-E.) (3. GF 
E E8 Q j= 1 J J 
‘ 
mn = 0 er, 2a. ; 
where 
De Z Pn vue 2.4, oe 
R. = IC, n + Cc, (n“-1)" + C, (en -l) Jw. 
2 
=| eo lean + een cap ) - C_n (en -l)v" (3. 25 
4 E j 3 J 
2 2 2 a 
-{C,en (2-en - €@°~] (w, 9 
4 tee 
and 
1 if ne=e6 
6 9 cd 
f: 0 ifu1 = 0 


The number of independent variables has now been re- 
duced from three (&, 8, and T) to two (& and T). However, 
now there are three equations to be solved for each value 


of n considered. 


Boundary Conditions 


The boundary conditions must also be expressed in 


terms of the Fourier coefficients of the displacements. 
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DMobyolpeOorted.. cubstitution of (3.2b)-(3.3a) into 
(2.24)-(2.27) yields the boundary conditions for a simply 


supported shell 


V = Q (3.12a) 
w= l= (0 (3-125) 
n 
= Q BAZ 
Ts ( c) 
ne = 0 at € = 0,1 


With (3.5a), @8.12a), an@ (3.12b) the last of the above 


Gonaditions reduces to 


nN Nn nN 
lure + Cyy + C 


n n = oe 
tx t V (Cram + C )} = 0 cece .0 51 


did ale 


From equations (3.12b), (3.4b), and (3.4e), it follows that 


n n 
Crom and Crop are zero at the ends of the shell. Hence, 


the final boundary condition becomes 


nN 


lur; ate ee one } = 0 at € = 0,1 ee kd) 


n 
ox 


Clamped. The boundary conditions for the clamped 
shell are obtained by substitution of (3.2) into (2.28)- 


(2.31). The result is 


vu = 0 (3.13a) 
vy) e=m0 (3.13b) 
a — (3.13c) 
wre = 0 at — = 0,1 (3.13d) 
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Response Symmetric to the Midspan. Substitution of 


(3.2) "ameo 22) = (2285) one es 


u = 9 (3.14a) 
ar = 0 (3.14b) 
Were = 0 (3.145 
Wreeg =a aul =goe 2 (3.14d) 


TTI MODIFIED FINI DItPERENCES 


In this section a brief discription of the differenc- 
ing procedure used to approximate the axial derivatives of 
(3.6)-(3.14) is presented. These finite differences have 
been called "modified" differences [31, 321 on “halt 
station" differences [33]. 

As an alternative to the finite differences, the de= 
pendent variables could have been expanded in Fourier series 
in the axial coordinate. T. Wah and W. C. L. Hu [34] made 
use of expansions of this type in their linear vibration 
analysis of simply supported, ring-stiffened cylinders. 

The equations of motion of the shell and the equations of 
motion of the rings were satisfied simultaneously by enforc- 
ing the continuity of displacements, forces, and moments 
between the shell and the rings. However, the present 
problem is more complex than their problem due to nonlin- 
earities in the equations of motion of the shell and due to 
the clamped boundary conditions. The use of the finite 


differences to approximate the axial derivatives appears 
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to facilitate the handling of both the boundary conditions 
and the effects of the discrete ring stiffeners. The use 
of series expansions in the axial coordinate was not in- 
vestigated in detail; consequently, the advantages and 
disadvantages are not clear at the present time. 

Modified differences are applicable to systems of 
equations in which only even or odd ordered derivatives 
of a given dependent variable appear in any one equation. 
As an example of equations that fall into this category, 
consider (3.6)-(3.10) with ne replaced by (3.5b) and 
without the nonlinear terms and the terms associated with 


mt ee See er SS Se 


the rings, that is 


n 1+v n 1-v 2 n n 


Uree + a ee ia 8 SEU, = 0 (3.15a) 

TY ==sO ow kee 
v = 

-n°yn — 2 ure + SY vi, tw = 0 (3.15b) 
Mos 2 

1 2 n Zee acen 

=a «|= (CW, - 2n w, Sew) (See) 

bz 6EEE Es 

= 6 p- ee a en vu, 

On 


Considering the spatial derivatives, note that only 
even derivatives of uw appear in equation (3.15a); only odd 
derivatives of u Gpecar time(s. 6b) and (3.15c). Similarly, 

n 


equation (3.15a) contains only odd derivatives of v" and w : 


and (3.15b) and (3.15c) contain only even derivatives of 


Or/ 


v and Wo For systems of equations having this property, 

a staggered finite difference mesh can be used advan- 
tageously. A staggered mesh, as illustrated in Figure 3.1, 
is one where half stations are interspaced between the 

whole stations, and the distance between consecutive whole 
or half stations is denoted by d. In the modified 
difference method, some of the dependent variables are 
defined only at whole stations and the remainder are defined 
only at half stations. For reasons which are not yet 
obvious, u" will be defined at half stations; while v and 


w will be defined at whole stations as shown in Figure 2. 


i 1/2 2 3/2, Kk? Loewe k k+1/2 k+l M-1/2 M 


a ; 
oe 


uy u. 
k-% k+35 
Ae cn ae 
oe k k+1 
va Te, ie 
k-1 k ae 
Figure 3.1 


STAGGERED FINITE DIFFERENCE MESH 


K. P. Chuang and A. S. Veletsos [31] were apparently 
the first to employ the modified differencing method in 
analyzing circular cylindrical shells» They proposecdmene 


fundamental modified central difference operator 
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EO haga. 


(3.16a) 


for derivatives with respect to the independent variable x. 


Consistent applications of the 


operator yield a set of 


difference approximations for the higher order derivatives 


of the form 


i! 
f, — eee - 2f, +f 3. LGb 
a lS aia et fy) ( ) 
Ll 
(4 ), = -—= (f, 3 - 3£ + 3f — —, (3). Hi@ie ) 
XXX k AS kts k+h k-% k- 
es = = 

(Er yea = aa (Th 42 40d - ae at zs f9) 

(3.16d) 


Returning mow to (3.15), 
obtained by expressing (3.15a) 
and (3.15c) at whole stations. 


equations are 


n 
u 


i n 
= + 
a2 k+3 eee ee 
ae ee en 7 uF 
2 k+% Va 
2 a ety i a 
n Vie 5 n 4 (arenas 
+ nw, = 0 
a n n 
1D a 4 (Wo Awe + Ow) 
+w ) + niw?} = 6 p - WU 
-] k On k k 
6 


the difference equations are 
at half stations and (3.15b) 


The resulting difference 


a sl cae saat n 
n — n = = 
(Wrod Wy.) 0 n Orly ae a « 
2 Ll-v 7 mn 
pe a — (Vi4y 2vy,. tv) 4) 
Ne = cchyeeeZe, (3.17b) 
As ie ee Soa 
k-l k-2 qe k+l k 
Cc “n ” n i a n = n 
in (OP IL (3.l/c) 
9 


Observe that in the development of these three equations, 


N at whole stations. Like- 


it was not necessary to define u 
ae ae 

wise, w and v. did not require definition at half stations. 

These circumstances prevail for two reasons: (1) the deriv- 


atives of any one variable are all odd or even ordered 


within a given equation, and (2) the equations are linear. 


Superiority of Modified Differences over Conventional 


Differences 

Chuang and Veletsos [31] showed that when the con- 
ventional central differences are used for all the deriv- 
atives, the difference equations obtained by discretization 
of the governing differential equations do not agree with 
those obtained by minimizing the finite difference form 
of the potential energy. On the other hand, when modified 
differences are used, the difference equations obtained by 
both procedures are identical. Consequently, the conven- 
tional and modified differencing of the differential equa- 
tions gives rise to two different numerical models, a 
conventional difference model that approximates only the 
differential equations, and a modified difference model 
that approximately satisfies both the differential equa- 
tions and the condition of minimum potential energy. In 
Appendix A, this inconsistency between the two schemes is 
shown to prevail even after the displacements have been 


expanded in Fourier series in the circumferential coordinate. 
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Chuang and Veletsos [31] compared the accuracy of 
the modified and conventional differences when used in 
the static analysis of cylindrical shell roofs. They 
found that the modified differences were far more accurate 
than the conventional differences for the same mesh 
Spacing. Ball [32] used both the conventional and modified 
differences to compute the natural frequencies of circular 
rings and compared the numerically computed frequencies 
with the exact frequencies. The frequencies obtained 
using modified differences were considerably more accurate 
than those obtained with conventional differences for the 
same mesh spacing. N. J. Cyrus and R. E. Fulton [33] 
developed analytical representations of the errors asso- 
Cliated with the modified and conventional differences 
when applied to beam and string problems. For a fixed 
number of stations, the errors in calculated deflections 
and curvatures resulting from the use of modified differ- 
ences were, in general, smaller than those obtained by 


using conventional differences. 


Adaptation to the Nonlinear Problem 
The discussion of the modified differences has avoided 


applications to nonlinear problems until now. 


Finite Difference Grid and Notation. In adapting 
the modified differences to the nonlinear equations, (3.6)- 
(3.14), the pattern used for the linear equations is 


followed. Difference equations are written which make use 


ae 


of a staggered node system, and the displacements are 
defined at either half or whole stations. A convention 

is adopted whereby both half and whole station quantities 
are given integer subscripts. To distinguish whole station 
quantities from half station quantities, an upper case 
letter will be used for the whole station subscript, and 

a lower case letter, for the half station quantities. 

As shown in Figure 3.2, the cr half station node is located 
between whole station nodes K and K+tl. The distance between 
adjacent whole station nodes, as well as the distance be- 
tween adjacent half station nodes, is equal to d. Let the 


radial and circumferential displacements w and v” be 


defined at the whole stations; and, the axial displacement 


n yh yh 
a YK VKa1 K+2 M 
ee we ee an wee 

il K teas KZ M 

Figure 3.2 


SPATIAL DISCRETIZATION USING MODIFIED DIFFERENCES 
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nN 


u , at half stations. 


operator of the modified differences becomes 


C ( aa = 


Qu] 


rit) 


k al 


and at half stations it becomes 


C ( rely = 


At whole stations the difference 


Here, and in sequel, k = K within a difference expression 


or equation. 


Difference Expressions LOG em DeR Evatt. Ves. 


; i) i 
three displacements u,v , and w 


nN 


Since the 


are not defined at every 


node point, the derivatives can be evaluated at either 


Wale statlLons Or 


derivatives that 


(Ving). = ‘ 
(Wee) = 5 
eek = 7 
(We eee) 


whole stations, but not at both. 


The 


can be evaluated at the half stations are 


il n 
ee (w 
Ae K+2 


n n m1 


(3.18a) 


(2 . 16) 


(3.18c) 


(3.18d) 


The derivatives that can be evaluated at the whole stations 


axe 


Mm -i 
(Uredy - 4 


(U4 1 7 Uy) 
k 
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(3.18e) 


n n n n 
(W,- ie = faa = 2W, 1+ Wee) (3.16) 
© Com a2 K+1 K K-1 
eS eee Fp Ss) (3.18g) 
Pee ae K+1 K Kaul ; 
n ee: n _ i nae n n 
(Wrrree) x = aa (Weis AW + Ow, AW 4 + Win 
(3.18h) 
Evaluation of Nonlinear Terms Mme nenlinearicrves 


in the axial and circumferential equations (3.6) and (3.7) 
are embodied in the functions Al” and Oe given by (3.8) 
and (3.9); in the radial equation (3.10) the terms Np Nps 
Ne@: and ng are nonlinear. The computation of these non- 
linear quantities in terms of the discrete displacements 
presents difficulties not encountered with the linear 
equations. 

To demonstrate the difficulty, assume that Oe must 
be computed at the whole station K. From (3.9), note that 
Cy? 
at station K in order to get Me Now, Co is defined by 


along with several other quantities, must be computed 


equation (3.4a), and is computed From 


ies © a i+n i-n 
C=, ae, (ne ve w! ! 
Oo. 726 E E E 
where 
Liem = 0 
&N — 
0O1if n = Q 


See Appendix B of Reference [30]. 


74 


1 aioe se N) 


2 eae = Th 


Since Co must be computed at node K, Wee must be known 
at whole station K for each value of n. As evidenced by 
(3.18b), the first derivative of w' can be expressed only 
at the half stations. Consequently, an interpolation of 
some form must be utilized in order to define Wre ata 
whole station in terms of w’ at nearby whole stations. 
The same situation exists for several of the other non- 
linear terms. The interpolations used to define the other 
derivatives falling in this category are developed in 
Appendix B. One example of these interpolations is the 

? 


expression for w at whole station K 


g 
n ee n _ 2 _ n we ue 


Mecation of Rings 


In the interest of simplicity, the assumption was 
made that rings could be located at whole stations, but not 
at half stations. As an aid in expressing the ring- 


associated terms, let the set H and the symbol See be 
i 


defined as follows: 


He = ae Sieh teinet oe = K, where K is the whole 
station location of the cp ring; j = 1,2,...,N} 


@3.. 99! ) 
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| Liga Hx. 


, ~ J (3.19b) 
a oO if Kk ¥ H, 


IV. DIFEERENCESD LEEERENT PAR SEOUAT IONS 


Hquaeroris “Or Motion 


In this section the difference-differential equations 
which constitute the numerical model are obtained by re- 
placing the €-derivatives in the governing equations by 
the modified differences of (3.18). 

When the response iS not symmetric with respect to the 
midspan, the finite difference grid is arranged so that the 
whole stations K = 1 and K = M are located at & = 0 and 
=", respectively. The total number of whole stations 
is thus equal to M. The half stations k = 0 and k =m are 
imaginary points located a half mesh space from the ends of 
the shell. The whole stations K = 0 and K = M+l are imag- 
inary points located a full mesh space from the ends of 
the shell. Thus there are M whole stations and m-l half 
stations on the shell. 

When the motion of the shell is Symmetric with respect 
to the midspan, the finite difference mesh is arranged so 
that there is a whole station at & = 0 and a half station 
at € = 1/2. The whole station at §€ = 0 is assigned the 
index K = 1, the half station at & = 1/2 1S assigned the 


index k = m, and m = M is now the number of half stations 


“See Pogue errs ac. 
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as well as the number of whole stations located on the 
half of the shell between € = 0 and € = 1/2. Thus, the 
whole Station K = M is located at € = *tv= Xd. 

By expressing (3.6) at the half stations, and (3.7) 
and (3.10) at whole stations, the difference-differential 


equations that govern the response can be obtained in the 


Beene 
Ny 
n uy e 22 ,1-v 7 2.22 2 
oe ene ———} ve {2 tdadn S =" n (Cn 
n sis n n 2 Zn 

+ C.) Sere Uy, dn ( ) eee u,_) = a Al, (3.20) 

N 

7 Neen? -aue eae ie 

=i “ 4 ii 5 Ki. Ten we 

) = J 

K=k=1,...,M-l1 (or M); n=0,1, 

3 


The rings were assumed to be located only at whole 
stations. However, in equation (3.20), the term 


r 
- 2 g2e2 (C ac + 3C)) UG 

is multiplied by u?, a displacement which is defined at a 
half station eos one-half of a mesh space from the ring 
at station K = H. Since the axial displacements are gener- 
ally small cated to the radial displacements, this 
inconsistency is probably not significant, especially if the 
grid points are close together. Alternatively, the above 
term could be replaced by 


N 


- & Aan (C n~ + C_)(6 + 6 ) 
= 
® ] ® +- H. 


Then half of the effect of a ring is felt at station k if 
there 1S a ring at station K or Ktl. 


7] 


it) n l+v n BPs a n 
(Foxe. 7am ere 2 tte int n 38 eH) Vi 
] Ey 1-v 2 
ss n ae n — n 
+ dn ( , ) Uy +: 5 a d 2, 
N 
a G32 19 
2 Z n 
a dA C, (en —-l) OTH. on 
eel J 
K=k=1,...,M=-1 (or M) ; n=O, > l= 
Tia oh n n n 
dj Wy an P, + dow, + dW q Wel + - Wy 
n n n 
+ dw te Ge teal) + n + 
4 K-1 3° K-2 O Ex a 
N 
ee Rs (32op 
8x y=. KH J 
K=1, ce, N=Opee oe 
Here Gyre sd, are Known quantities which depend onon, c, 
n 
a, and d. The terms es Dene ie ; oi : no eand are 
co bo 


calculated as explained in Appendix B. In (3.20) and 
(3.21), the upper limit M of the values assumed by K applies 
only when the response iS symmetric with respect to the 
midspan. The ring term Re, obtained from (3.11) and (3.18), 


is given by 


RE = d wey + dow + dew) + dg(up - uy) + d vy 
Kees J=1,2,...,N) (322559) 
Here Gass++1d, are also known quantities that depend onn, 
d, and the ring coefficients Cyreee Cee and e. The expres-= 


Sions for dj reese rd are given in Appendix C. 
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Boundary Conditions 

Casting the boundary conditions in finite difference 
form is conceptually simple, but the process involves 
development of a few interpolations. The interested 
reader will find the development of the interpolations in 


Appendix B. 


Stipe lyY Supported. The bourdary conditions (3.12) 
allow an axisymmetric rigid-body translation along the 
E-coordinate. Consequently, when only the conditions of 
(3.12) are observed, the numerical solution is not unique. 
In order to obtain a unigue solution, the axial displace- 
ment of the axisymmetric mode must be set to zero at some 
arbitrary axial station. To facilitate a numerical 
solution, the node at &€ = 1 where K = M was chosen. Since 
the axial displacement is not defined at the whole stations, 
interpolation is required to express this condition. The 
interpolation is developed in Appendix B. Substituting 


the appropriate equations from (3.18) into (3.12) at K =] 


yields 
wy = 0 (3.24a) 
w} - 2w, + wy = 0 (3.24b) 
vi = 0 (Ge2ae) 
lie mp at = 
uf - up t+ A(Cyy, + Crxx,) 0 (Groce) 
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and at K = M yields 
Wie 0 (3.25a) 
We ~ 2Wy + Wy = 0 (3.251) 
a = 0 ( 3ia,2 Ses) 
Un -~ Un? + dCKa . CTxx,) a0 (3.254) 


In order to tie down the axial displacement of the axi- 


symmetric mode,the following expression is used at the end 


4 
Saal 
0 0 iPad 7 
u, + oT = A v= = 0 (3.26) 
where the half station k = m is located at —& = 1 - sd. 


Clamped. The finite difference formulation of boun- 
Gary conditions s(@813b) andy (3F13c) 1s strargqntrorwara 
requiring only the substitution of the appropriate expres- 
sions from (3.18). However, the discretization of (3 2igay 
and (3.13d) requires the utilization of interpolation form- 
4 


ulas. The difference expressions for the clamped boundary 


conditions are 


n n ens 

Uo + 2u, wel Us = 0 (3.27a) 
V7 = 0 (3.27b) 
wi = 0 (3.27c) 


See Appendix B. 
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Wo - 9w, + wm O Nek = (aout 


3 
and 
De Oi ee (ales) 
Un F “Unig 7 3 Ym-2 = .28a 
ieee 
Vu = Q (3.28b) 
Way = (3.28c) 
n n n = Ee: 
Wy 7 IW + Wut1 = 0 at K=M (3.28d) 


Response Symmetric to the Midspan. Substitution of 


the appropriate expressions from (3.18) into (3.14) leads 


EO 


ue 0 (3.29a) 

Vin > Ved =_0 (3.29b) 
a ee 2 

Wy Weed = (3.29c) 
nN mm al 
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CHAPTER IV 
SOLUTION OF THE NUMERICAL MODEL 
I. SOLUTION PROCEDURE 


There are various methods that could be used to solve 
the differential-difference equations that comprise the 
numerical model. Several of them are examined in this 
chapter. The algorithm that appears to be the most promis- 
ing 1S a combination of the Gauss elimination method and 
the iterative Newmark beta-method. The solution procedure 
1S carried out in the following manner. Suppose that Weer 


-n vor nm n : 
Wer Wer ule and Vy, are known at time Ty for all values 


of K, k, and n. The solution is sought at time Toe a small 


increment later than tT): In the first iteration for the 
: een ; “on 

Selution at Tos ee at tT, is taken equal to Wr at tT): The 
ae rel 

quantities We and W, at T, are then calculated for all K 


and n uSing the Newmark beta-method. Knowing WK at T the 


5! 
terms es and Noy which depend only on Wier are computed; 
and oe and Ve are obtained by solving equations (3.20) and 
(3.21) with the Gauss elimination method. Next, having 
computed estimates for Wye Wan Be and Vig at time Tos the 
acceleration vy at T, can be calculated from equation (3.22) 
for all K and n. The accelerations of the two consecutive 


iterations are compared to check for convergence. The 


iteration is continued until convergence is achieved. 
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This procedure is repeated for many time steps to compute 


the displacements as a function of time. 
Ti. SOLUTION OF THE AXIAL AND CIRCUMFERENTIAL EQUATIONS 


Applicable methods 

The axial and circumferential ordinary differential 
eqjuacions, (3.60) and (3.7), Censtitute a linear boundary 
value problem in the two unknowns ee and 7 te) when Mele 
ee and w’ are known. In the development of the numerical 
model, these equations were replaced by a system of simul- 
taneous algebraic equations (3.20) and (3.21). If Aly 


n ; 
A2 and Wy are known, these equations can be solved in- 


n 
K! 
dependently of the radial equations (3.22). Either a 

direct elimination procedure or an iterative procedure is 
normally employed to solve systems of algebraic equations 

of this type [35, 36]. An example of the direct procedure 
is the Gauss elimination method. This method has been used 
by several investigators to solve static shell problems 

[23, 30, 37]. The Gauss-Siedel method is a popular itera- 
tive procedure. 

The Gauss elimination method and the Gauss-Siedel 
iteration procedure were compared in an attempt to determine 
the best of the two methods for the problem under considera- 
tion. Computer routines using both methods were written 
to solve (3.20) and (3.21). Solutions to some simple 
linear problems were generated by both routines, and the 


computed solutions were compared with the exact solutions 
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to the differential equations. The Gauss elimination pro- 
cedure required Significantly less time on the computer 

than did the iterative procedure for the same degree of 
accuracy. When both methods used the same execution time, 
the elimination method provided a considerably more accurate 
solution. Consequently the Gauss elimination method was 


selected for the algorithm. 


Gauss Elimination 
To facilitate the development of the Gauss elimination 


procedure, the axial and circumferential equations are 


converted to matrix equations. First, let Aly and 2 be 
defined by 

Ny 
Al, = qa? ere 4- - KH, an® [C, (en*=1) - Cyl (weiy - wy) 

Ne (4 .la) 
A2n = a? M2, + pw Sa, an C, (en*-1) Wy (4 .1b) 


Equations (3.20) and (3.21) can now be written in the form 


nf Nn 
“ae 7, = 


FT iA 
K K-1 Jy 


A Z a BD 


rN 
peal K* (4.2) 


where n = 0,l,...; and K = 1,2,...,M-1 or M, dependingses 
whether the response is unsymmetric or symmetric with 
respect to the midspan. For the unsymmetric response, M is 
the number of whole stations between € = 0 and € = 1. If 
the response is symmetric with respect to the midspan, M 


is the number of whole stations located between € = 0 and 
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E=1/2. The matrices ae Byr and Gu are 2x2 matrices whose 


elements are given in Appendix C, and Ze and ox are defined 


by 

yy? 

a ? (4.3) 

a = _ : 
“kK 
a 

4 k 

oa (4.4) 
A nN 
= 2k 


For a fixed value of n, (4.2) represents a system of 2M-2 
(or 2M) equations in 2M+2 (or 2M+4) unknowns for unsymmetric 
(or symmetric) response as shown in Figure 4.1. The bound- 
ary (or symmetry) conditions provide the remaining neces- 


Sary equations 


Recursion Relations. The Gauss elimination recursion 


relation for equations having the form of (4.2) is [30, 37] 


n nN 


= PK 2K4]1 


n n 
Ze + X; (4.5a) 
where ae 1s a 2x2 square matrix, and 7 is a 2xl vector. 


Thus, the recursion relation for 2a is 


Zeiy = 7 Pei 2 + XRH-1 (4.5b) 
Substitmrien ter (2.58 )"4into (4.2)* results’ in 
n _ n Ween a een n al 18 =A us non 
Zz, = ~(Bp - C Pyy) A 2, +(BK - © Px-1]) Igx-C  xXy_4] 
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€=0 AXIAL EQUILIBRIUM €=] 
. SATISFIED AT 
2] 3 HALF STATION NUMBE m-2 m-l | m 








O 
0 = -X— -O —x-——-O—  —O— 
0 CIRCUMFERENTIAL 
EQUILIBRIUM 
SATISFIED AT 
ws 4 i 3 WHOLE STATION NUMBER) ME OM 


EQUATION (4.2) 
(a )UNSYMMETRIC RESPONSE 





€=0 € 2 AXIAL EQUILIBRIUM 
a ee SATISFIED AT 
HALF STATION NUMBER 
O=x— x—o—x-f X——=Q——_X- Qe —X 
5 M+l ( CIRCUMFERENTIAL 
EQUILIBRIUM 








SATISFIED AT 


i Cy <a Ro WHOLE STATION NUMBER 
~ 2 


M-| M“®-INDEX K OF EQUATION (4.2) 


(b) SYMMETRIC RESPONSE 


FIGURE 4. 


RELATIONSHIP BETWEEN MESH NODES 
AND 
INDEX K OF EQUATION (4.2) 
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Comparing the above equation with (4.5a) reveals that 


nn cal 


nN nN Nn 
= seein A 4.6 
Py IB Gs ot ( ) 
and 
De tet - pty tg - co? 0) 4.7) 
en ey on Xe-1 (4. 


In solving (4.2), the following steps are performed: 
(1) the Px matrices are computed by recursion relation 


(4.6); (2) the x2 


K vectors are calculated from (4.7); and 


(3) the solution vectors Zi are computed in a backward 
sweep by recursion relation (4.5a). Each of the steps 


will now be examined in more detail. 


Computation of Py Matrices. The general matrix Pe 
depends on ae Hence, Py must be determined first. This 
matrix is determined from application of the boundary con- 
ditions. The boundary conditions which uniquely determine 
the Py matrix are given by equations (3.24c) and (3.24d) 
when the shell is simply supported and by (3.27a) and 
(3.27b) when the shell is clamped. From the two appropriate 
boundary equations and the axial equilibrium equation at 


K = 1 an equation is obtained of the form 
— pn : 


The elements of the matrix Py for both the simply supported 
and clamped shell are listed in Appendix C. By using the 


recursion formula (4.6), the remaining Py matrices are 
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computed for all of the participating modes and for 
K = 2, 3, «..-,M-2. When the response is symmetric with 
respect to the midspan, K = 2, 3, ...,M-l. These P 
matrices do not change with either the load or the 


response. 


Computation of Xi Vectors. When the shell is simply 


supported the vector xy of (4.8) is given ce 


n n n 
Aly + wae + Te, 
he | ie ac 2 1l-v 
x, = [ n eeu 
F 0 (4.9) 
and when the shell is clamped, by 
hee 
n : 2 
2 ro 
eit /(4 + don p)) (4.10) 
0 


Since in both cases x4 contains nonlinear terms in the 
radial deflection, x4 is not constant, but instead, varies 
with the load. The vector cae defined by (4.4), also 
changes with the load. Thus, the vectors Xi and ae must 
be recomputed many times as a solution is generated; x4 


n 
from (4.9) or (4.10), and the other Xp vectors from (4.7). 


‘the fact that the boundary condition on v, is vi = 0 
makes it unnecessary to introduce vf into the solution and 
to consider the circumferential equilibrium equation at 
point K = l. | 
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Back Substitution. The solution vectors cae are ob- 
tained through back substitution, that is by sequential 
solution of equation (4.5a). However, relation (4.5a) 
requires a "Starting" z-vector. 

The boundary conditions at §€ = 1 provide the infor- 
mation necessary to compute this vector when the response 
is unsymmetric with respect to the midspan of the cylinder. 
The boundary conditions are given by (3.25c) and (3.25d) 
when the shell is simply Supported and n # 0. When the 
shell is clamped, or when the shell is simply supported 
and n = 0, the corresponding boundary conditions are ex- 
pressed by (3.28a) and (2.s210 oe By substituting the 
boundary conditions into (4.2) with K = M-1 and by making 
use of (4.5a) at K = M-2, the following result is obtained 
for the solution at K = M-l. 


= ie oer ert 


n "ae ne ani 


O" Xy-9] (4.11) 


me 
Ms 
n 


nand £" are given in Appendix C for the simply 


where H , O 
Supported and clamped boundary conditions. 


When the response is symmetric with respect to the 


n 


M 1s determined by using the 


midspan, the starting vector 2z 
symmetry conditions given by (3.29a) and (3.29b). If 
(3.29a) and (3.29b) are substituted into the circumferential 


equilibrium equation at K = M and the result combined with 


(4.5a) at K = M-1,° the result is 


“Equations (3.26) and (3.28a) are identical when n=0. 


3the axial equilibrium equation is identically satis- 
fied at the midspan. 
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n n Ty rT n 
zy = (1 +R Di] Ree eee Sian (4a 


where I is the identity matrix and R and s. are given in 
Appendix C. 

Once equations (4.11) or (4.12) have been evaluated 
for each n, the remaining z-vectors are computed in sequence 


using (4.5a). 
III. TIME INTEGRATION 


Four basic methods were considered for carrying out 
the time integration of the radial equations of motion. 
They were the Runge-Kutta methods, the explicit methods, 
the implicit methods, and the Newmark beta-method. 

The Runge-Kutta method [38] has been used by several 
investigators to compute the nonlinear dynamic response 
of cylindrical shells [ll, 20, 22]. However, they included 
only a few modes in the response and did not employ finite 
differences along the axial coordinate. Consequently, 
their systems of equations were much smaller than the system 
to be solved here. For the large system of equations of the 
present analysis, the Runge-Kutta method would require 
prohibitively large amounts of computer storage space and 
excessive computing time. | 

The explicit methods [39, 40] all have the common 
feature that the unknowns can be calculated explicitly 
without solving a large system of simultaneous equations 


each time the solution is advanced one time increment. 
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Nonlinear equations can be solved without iteration since 
the nonlinear terms contain only products of previously 
computed quantities. However, the explicit methods are 
numerically unstable if the time increment is too large. 
While the largest permissible time increment can be pre- 
dicted theoretically for linear equations, the permissible 
time increment for nonlinear equations has not been estab- 
lished as yet. Since a numerical instability might easily 
be mistaken for a physically unstable response, the ex- 
plicit methods were ruled out. 

The implicit methods [39, 40, 41] all utilize time 
averaged finite difference approximations of the spatial 
derivatives. As a result, the difference equation at any 
spatial location contains several unknowns. The advance- 
ment of the solution by one time increment is accomplished 
by solving a system of Simultaneous algebraic equations. 
If the differential equations are nonlinear in the spatial 
derivatives of the dependent variables, as in the case for 
the equations under study, the simultaneous algebraic 
equations are usually nonlinear. These nonlinear equations 
must be solved by an iterative procedure that solves a 
sequence of linear equations. The implicit methods have 
the advantage of being numerically stable when applied to 
linear problems. However, since an iterative procedure 
must be used in the implicit methods to solve the nonlinear 
equations, the decision was made to use an iterative inte- 


gration procedure. 


91 


The iterative integration procedure selected as the 
most suitable for carrying out the time integration is the 
Newmark beta-method. References [42] and [43] contain a 
detailed description of the method and a theoretical deter- 
mination of the stability and convergence limits. A typical 
application of the method, as well as references to several 
other investigations that employed the method may be found 


in Reference [44]. 


Inte mo bdo nro ce alume 


The numerical integration procedure begins with the 


n 
assumption that the displacements, ae Ve, and w.; the 


K! 
- n . 

velocities, wy; and the accelerations, We are known at 

time T,- The corresponding physical quantities may be 


computed at time T, = tT, + At, where AT is the time incre- 


2 
ment, by repetition of the following steps: 

1. Assume that WK (TS) = WK (T)) Fon thea tia 
iteration. For subsequent iterations, assume 
that We (To) = Wy (T.) where wy is computed 
in step 6. 


2. Calculate the radial displacements and velocities 


at time tT. by the relations 


2 
on -n gal wel 
Wy (t,) = i + alw,(T,) + We (t,) | At (4.13) 
Wy (tT) = WK (T)) We (T,) At + (%-8) Wig (7) (AT) * (4.14) 


nn Z 
+ Bw, (tT, ) (AT) 


22 


The significance of the parameter 8 will be 
discussed in a later paragraph. 

3. Evaluate the nonlinear terms Al; and WW from 
(APT ai. 


n a 
4. Compute u, and VK by Gauss elimination. 


°K 


6. Compute the radial accelerations from equation 


n n 
5. Evaluate n, , Nes Neg: and Ng 
K K K 


(3.22). Let wl (Ts) denote the values computed 
here to distinguish them from Wig (T) of stepL 

7. Compare wly (Tt) with tix (15). If their relative 
difference is larger than a prescribed value e, 
the convergence criterion, repeat steps 1-7. 
If the difference is smaller than ¢e_, the solu- 
tion is said to have converged. 

This procedure can be repeated as often as desired 

in order to determine the response of the shell as a func- 


iaL@n sof time. 


The Parameter Beta. Equations (4.13) and (4.14) are 
the quadrature relations of the Newmark beta-method. 
Newmark [43] has shown that the choice of a particular 
value of 8B is related to the shape of the acceleration- 
time curve during the time increment At. Acceleration- 
time curves associated with g = 1/4, wien and 1/8 are shown 
in Figure 4.2. For this study, 8 was taken equal to 1/6 
corresponding to a linear variation of the acceleration 


during the time increment. 


a3 


wl 





Figure 4.2 


ACCELERATION VS. TIME FOR 8 = 1/4, 1/6, AND 1/8 


The stability and convergence limits of the Newmark 
beta-method are given by critical values of the ratio of 
the time increment ATto the shortest natural period T of 
the finite difference model. If the ratio AT/T is larger 
than the stability or convergence limit, the procedure is 
unstable or does not converge. The stability limits, 
derived by Newmark [43] for a linear oscillator, are given 
in Table I for several values of 8. When 8 is larger than 
one-eighth, the convergence limit is more stringent than 
the stability limit. 

The problem under study is governed by ee of gnene 
linear equations, and the shortest natural period is not 
easily obtained since the natural frequencies of a non- 


linear system vary with the amplitude of the motion. 
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Nevertheless, the shortest natural frequency of an un- 
stiffened shell, as predicted by the linear theory, will 


serve as a useful guide. 


TABLE I 


STABILITY AND CONVERGENCE LIMITS 


B 0 Paez 1/8 1/6 1/4 


Stability Limit, 
At/T = 0.318 0.389 0.450 0.551 infinite 


Convergence Limit, 
ATtT/T = iin Ono ome wr Od 5:00) 0x38 OMe ars 


For the Fourier modes of interest, the linear equa- 
tions reveal that the axisymmetric mode (n=0) has the 
shortest ae If the shell is simply supported, 

; = ne@= O and ny = Pe aa Hence, from (3.10), the 


linear differential equation for the axisymmetric mode is 


Nn 


seen to be 


0 


0 2 0 : 
ee fede > ) wove (4.15) 


OQ 


12 WYEEEE 


Approximation of the spatial derivatives by central differ- 


ences leads to 


- 0 0 0 0 0 
& aa hs 
=e (Weao 4Wyiy + Wy 4wy_1 + Weg) (4.16) 


{3 
Sle 


P(l=v- yw, = -w, 


4mn Chapters vt will be shown that this is true. 


25 


Assume a solution €e, (4.15) eel enero 


. ; qné 

w = Sin ae Sin ae. (4.17) 
where Og is the natural frequency of the axisymmetric 
mode with wavelength A If the shell is divided into 
M-l intervals of length d, it follows that 

E Sak 1) ee 

1 = (M-1) a 
Thus equation (4.17) may be written 

wo = sin (8, 1) sin (WK) (4.18) 

Og M-1 





The natural frequency si is now the frequency associated 
q 





with the difference-differential equations. Substituting 
(4.18) into (4.16) and performing several algebraic and 


trigonometric manipulations yield 


_ dean” 2G 2 qt ae 
ave = oR ry. [2 cos a 8 cos (oD + 6] + (l-v ) 








Hhaesmea imum voles a occurs when ais = 1. Hence, the 


q 


highest frequency iS given by 
+ a (4.409) 


Consequently, the shortest natural period of the finite 
difference model is 


2 =5 
> 40 8 
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With T given by (4.20), At for convergence and 
stabm@lity is ~iven in Table I as a function of 8. Thus 
taking 8 equal to 1/6, 


4 
At < (0.389) (27) (> — + 1l-v ) (4. 2a) 
e a4 
insures both convergence and stability of the algorithm 
in computing a linearized solution. The assumption is 
made that the same reasoning can be applied to the non- 
linear equations; that is, if the calculation converges, 


it is also numerically stable. 


Convergence Rate and Computational Errors. The rate 
of convergence used by Newmark is defined as the error in 
the derived acceleration divided by the error in the assumed 
acceleration. Newmark has shown that for a linear oscilla- 
tor, the rate of Saye eas as well as the error in the 
amplitude and frequency of the computed motion, depends 
on the size of the time increment. These relationships 
are shown in Figure 4.3 for 8 equal to 1/6. Inspection 
of Figure 4.3 reveals that in order to keep the errors 
small, the time increment must be much shorter than the 
convergence limit permits. For example, if errors in 
amplitude and frequency larger than two percent can not be 
tolerated, then AT/T must be less than 0.12, whereas the 
convergence limit requires only that AT/T be less than 


OOo . 


oT] 


LO} - 1.0 





convergence 
limit 


2 
Ee 
o 





RATE OF CONVERGENCE -#2 
RELATIVE ERROR IN PERIOD~ #1 


FIGURE 4.3 
DEPENDENCY OF CONVERGENCE 


RATE AND COMPUTATIONAL 
ERRORS ON LENGTH OF TIME 
INCREMENT 
B = (46 (Ref. [43]) 


RELATIVE ERROR IN MAXIMUM RESPONSE TO AN INITIAL VELOCITY - #3 





OOIL-| .Ol OO 
O 10 .20 30 .40 
ATtT/T 
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Newmark [42] has suggested that when 8 is equal to 
1/6, a stable integration procedure for nonlinear equations 
can be obtained by requiring that the assumed acceleration 
and the derived acceleration agree to within one percent, 
or less, by the end of the fourth iteration. If conver- 
gence is not obtained by the fourth iteration, the time 
increment should be reduced. On the other hand, if the 
accelerations converge in one or two iterations, the time 
step should be increased to ensure optimum computing 


efficiency. 
V. COMPUTER PROGRAM 


A computer program was written for the IBM 360 Model 
67 computer to generate a solution to the governing 
difference-differential equations developed in Chapter II 
and Chapter III. Appendix D contains a description and 
complete documentation of the program, including instruc- 
tions for the preparation of input data cardsand other 
information necessary for the operation of the program. 
The program is dimensioned to accommodate twenty equally 
spaced axial stations (Ms 20) ane fifteen Fourier modes. 
It will handle up to three ring stiffeners located so as 
to divide the shell into segments of equal length. Both 
Simple and clamped supports can be handled. The basic 
output consists of the displacements and in-plane forces 


at specified time steps. The axial location, time, and 


og 


magnitude of the maximum circumferential stress at @ = on 


are also printed. 

The program adjusts the size of the time increment so 
that the number of iterations always falls within a range 
specified by the user. The number of iterations required 


to obtain a converged solution depends on the convergence 


criterion, €. Convergence is assumed to have occurred when 
Ti I 
coil esl) 
| (Wag) a (wi) | De 0. jlo eae 
E er 
- I+] i 
K = 0,. M 
“ND oN peoeer 
Max | (aR) |, | (WE) } 


whee f£ is the itenpambons numer 


Calculation of the Circumferential Stress at a Typical 
ocak rom 

A structure designed for repeated use may experience 
low cycle fatigue cracking and even total failure if the 
stresses repeatedly exceed the yield point stress of the 
Macerial. Therefore, one criterion, pom sa ssateelecemals 
that the stresses which develop during the response must 
never exceed the yield stress. For the problem at hand, 
only radial loads are considered; and the load is carried 
principally by circumferential stresses. Thus, the circum- 
ferential stress is a close approximation to the maximum 
principal stress. In the present program the circumferen- 


: 5 
tial stresses are computed at 0 = Cheat 


“At a fixed axial location, the circumferential stress is 
not necessarily largest at § = 0°7 but Since the assumed 
imperfection modes are in phase at 6 = 0%, it is a reason- 
able indication of the state of maximum stress in the shell. 


OO 


At a point on the midsurface of the shell, the cir- 


cumferential stress S%, is given by 





o, = 7 Ng (4.22) 


At either the inner or the outer surface the magnitude of 
the circumferential stress is larger than the stress at 
the midsurface, the additional stress being due to bending 
of the shell. The circumferential surface stress due to 


the circumferential curvature change, denoted Sop: is given 


by 


(4.23) 


[5 
ho} 
OQ 


Q 
Il 
I+ 
NS) 
a1) 
OQ 
cD 


NO 


OB 


The influence of the meridional curvature change on the 
circumferential bending stress is neglected. The magni- 


tude of the largest surface stress, denoted |co. |, is 


Qs 
given by 


Jon. = Roe Io gp! (4.24) 


Substituting (4.22) and (4.23) into (4.24) leads to 








lo (4D 


gs | 


On substituting (3.2c) and (3.3b) into’(4.25), the result 
is 


ii co 
ng cos né| + = Eal £ n“w" cos né| (4.26) 
0 n=0 





im g 


acO 


= Ea| = nw" | (4929) 
n=0 





Let n, be defined by the relationship 





9 
EB —_ 
fon | = n (4.28) 
Os g=0° 1-v2 0 
Then from (4.27) and (4.28)° 
en | etl =v— ele aaa (4.29) 
9 ae 5 e 
n=0 n=0 


At a given axial station, the value of Ny is thus propor- 


tional to the stress at 6 = 0°. Equation (4.29) can also 


be expressed in finite difference form as 


nN, = | x ng | + 5 a| = nowe | (4.30) 
K n=0 K n=0 
The maximum circumferential stress at 0=0°, Np , OCCURS 
MAX 
at some grid point K. Hence 
Np = ax |i, me Kee 23 7 ee, Ma eo mt (4. 38) 
MAX x 


For convenience, we will refer to this as the maximum 


circumferential stress. 


orn the computer program, the quantity ee in (43235 
was taken equal to unity. The error introduced is cer- 
tainly no greater than the errors involved with the other 
approximations. If the program is revised to compute the 
exact stresses, this factor should be reinserted. 
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Verification of Computer Program 


As indicated in Chapter I, there apparently are no 
previous solutions for the nonlinear dynamic response of 
a cylindrical shell that can be used to check all of the 
features of the present computer program. Furthermore, 
the program cannot be readily simplified or altered to 
correspond with the previous analyses. Consequently, only 
a linearized version of the program was verified through 
comparison with other solutions known to be correct. It 
might be argued that some of the new features of the pro- 
gram are of minor significance, and that the results of 
some of the earlier analyses could be used to test the pro- 
gram. However, if discrepancies were discovered, it would 
be difficult to decide whether to attribute them to an 
error in the program or to differences in the analyses. 

To check the capability of the program to compute a 
fully nonlinear, nonaxisymmetric solution, a typical 
problem was run through several time steps, and a complete 
sequential printout of the program variables was obtained. 
The printout was then checked against a hand calculated 
solution 

To date, those features of the program that compute 
the contributions from the rings have been checked ina 
cursory manner. One computer run was made for which the 
Shell was assumed to have ring stiffeners. The computed 


displacements appeared to be correct. The parts of the 


Ors 


program that compute the ring contributions are isolated 
and do not affect results obtained for unstiffened shells. 
The results presented in the next chapter also serve 
to verify the analysis and the program in that many of the 
results agree very favorably with those obtained by other 


investigators. 


Linear Test Cases 

Under an axisymmetric load, the response of the shell, 
as governed by linearized equations, is also axisymmetric. 
The correctness of the program in computing a nonaxisymmetric 
response was checked in conjunction with the verification 


of the nonlinear features of the program, 


Moving Step Load. The transient, axisymmetric re- 
sponse of a simply supported shell subjected to a moving 
step pressure discontinuity has been determined earlier 
by P. G. Bhuta [45] using a series solution to the linear 
Donnell equations. A subroutine was assembled to compute 
p(&, t) for this loading condition, and a numerical solu- 
tion for the axisymmetric response was computed using the 
program with all nonlinear terms set equal to zero. 
Twenty-one axial stations were used,’ and convergence to 
one per cent in a maximum of four and a minimum of three 


iterations was specified. A comparison of the radial 


the present program is dimensioned to handle a 
maximum of twenty axial points. 
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displacements ut the quarter-span station 1s shown in 


Figure 4.4. The constants for” the “example are also givem 
imeFigure 4.4. ‘The computed resilts ag®.:e fairly welll waren 
the series solution. However, by the end of the calcula- 


tllomnwawdetbectable error has®developed. 

The computed displacements are given in Figure 4.4 
Dino vwmeay tGmth tLinesstep. % Duringsthe time Ghat thes Read 
GLScOmMEINUELY wasm@em the Shell, hes time 2pemenenitgeemai ped 
small, and approximately one hundred and seventy time 
steps were taken to compute the first cycle of the motion. 
SHomie, berone the Giscemtinusty movedszofrLi ihe shell, the 


time increment increased, and approximately thirty time 


steps were required to compute the second cycle. The 
ratio At/T, where Ty i@ @eéiimed’ by equatiom (4.20), reached 
@ maximum vailue of .05. From Figure 4.3, it is seen that 


.he computational errors should be very smail for this 
Value. The numerical errors are apparently accumulative 
nd may reach signilticant magnitudes in lengthy compu- 
tattons. In subseyuent computations, the parameter € was 
eeuecqual to .00Meto reduce: the accumulative error. 
Uniform Harmonic Load. The axisymmetric transient 


response of a simply supported shell subjected to a loading 


given by 
Pia = Pp, cos QT 
was computed and compared to the series solution. The 


response predicted by the linear theory is 


Oe 


o Computed Points 
0.0010 — Ref, [45] 
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FIGURE 4.4 
QUARTER-SPAN REPONSE TO A MOVING 
STEP LOAD Ps 
(a= 18in, pie /2in, L=72in, Ay 20psi, E=30xl0 
p*7.332x10 Ib-sec in, vy =03, velocity = 72,000 in/sec) 
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oul 4p, 


——_—=——=— (cos 2T - cos 2 _T) sin (ams) 
Ceo Oper To ee ) q 1 


Thirty-one axial stations were used, ® and the convergence 
parameter was set equal to .001. The computed solution 
at every 10th time step and the series solution are shown 
in Figure 4.5. As with the previous test case, a small 
error developed after several cycles of the oscillation 
and was manifested primarily as a phase difference between 
the computed and exact responses. 

On the basis of the above tests, it was concluded 
that the algorithm provides an adequate representation of 
the axisymmetric transient response of a shell subjected 


to axisymmetric loads. 


Spuring the initial phases of the investigation, the 
program was dimensioned to accomodate 31 axial stations. 
The dimensions of the computer program were later reduced, 
and the program listed in Appendix D will handle a maximum 
of twenty axial stations. 


ON 


FIGURE 4.5 
RESPONSE OF MID-SPAN TO A 
UNIFORM HARMONIC PRESSURE 
(n/a = 1/36, L/a=4, Pe =.00218, 
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CHAPTER V 


NUMERICAL RESULTS 


Lt.) INTRODGE LIST 


The analysis, the numerical model, and the method 
of solution developed in the preceding chapters are appli- 
cable to a rather broad class of problems. A complete 
Study to determine the significance of the many features 
incorporated in the analysis and the computer program is 
beyond the scope of the present work. Hence, a problem 
that has been considered in two very recent investigations 
[19, 20] was selected for an in-depth study. Treatment 
of the problem by the algorithm of this dissertation ac- 
complished three things: (1) increased understanding of 
the solution to the problem was developed; (2) the feas- 
ibility and practicality of the new approach were 
established; and (3) deficiencies of the program were dis- 


closed. 


Statement of the Problem 

The problem is the determination of the dynamic 
response of an unstiffened, simply supported cylindrical 
Shell subjected to an exponentially decaying, uniform, 


radial pressure pulse. 


Shell Properties. The shell was assumed to be made 


from 6061-T6 aluminum and to have the following properties: 
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a ieee late) 


L/a = 2 
h/a = 1/100 
9 = 2.54 x 10°"lb - sec* - in’ 
E = HO psi 
Vv = = ie 
Op = 42 ksi 
c = 3.416 x 107? 
where ae 1s the yield point stress of the material. The 


nondimensional damping coefficient c corresponds to a 
value determined experimentally by Fung, Sechler, and 
Kaplan [46]. Yao [6] used the same coefficient in his 


adynamic Stability study. 


Initial Imperfection. The existence of a small, but 
finite, initial imperfection in each mode was postulated. 
The Fourier coefficients of the initial imperfection were 


assumed to be of the form 


Fal eae (7) Gel 


where wi’ is the amplitude of the imperfection in mode n. 
Except for one special case, the amplitudes of the imper- 
fection were all taken equal to 1/1000th of the shell 
thickness. The corresponding nondimensional amplitudes 


were thus 


ileal 


For the one special case, the imperfection amplitudes were 


taken from a table of random numbers. 


Load. The loading was assumed to be 


0 jeie ie << (0, 
P(x,t) = i (5 29 


Poe EG ioneese | 10, 
where Pe. is the peak pressure and ty is a time constant. 


The loading is expressed nondimensionally by 


0 fie ee) 


p(€,T) = (5.335 


Pace one) Pee 


For this study, a awameiie By integrating the load from 
t = -~m to t = +~, the impulse per unit area is found to 


be 


I=Pt (5. 


Ti. PREVIOUS INVESTIGATIONS OF 2HE Peeler 


Certain aspects of the problem outlined above have 
been the topic of two very recent articles. Anderson and 
Lindberg [19] investigated dynamic buckling under an 
exponentially decaying, uniform, radial pressure, given by 
(5.2). McIvor and Lovell [20] computed the dynamic response 
for a uniform radial impulse, a limiting case of the loading 


considered here and in Reference [19]. In the following 


ale 


two sections, these investigations are reviewed in some 
detail in order to lay the foundation for later comparisons. 
The simplified equations used by Anderson and Lindberg 
are examined in considerable detail as they illustrate 


many of the essential features of the problem. 


Anderson and Lindberg's Investigation 

Anderson and Lindberg's analysis identified the peak 
pressure and the total impulse as the important load para- 
meters for determining whether dynamic buckling will occur. 
When the total impulse was relatively small, they found 
that buckling occurred at very high pressures; in fact, so 
high that buckling took place only after the onset of 
plastic flow. A tangent modulus model was developed for 
predicting buckling loads in this regime. For loads of 
large total impulse, buckling occurred at much lower pres- 
sures and was governed by an elastic model. Both of these 
models were developed for the prediction of dynamic buckling 
loads and were not intended to predict the dynamic response 
under sub-critical loads. Their elastic model is of 
interest here and will be used as a standard of comparison. 

For their elastic model, they assumed the motion to 
consist of a uniform radial displacement plus an infinite 
number of nonaxisymmetric flexural modes. The flexural 
modes were required to satisfy simply supported boundary 


conditons, and the axial variation of both the flexural 


eS 


modes and the initial imperfections was taken to be a half 


Sine wave. These assumptions were satisfied by 


wo sin (HE) cos né 
1 i 


= 

i 

= 

+ 
lm 8 


for the radial displacement and 


co 
w= ¢ wi’ sin (Ty cos né@ 
for the initial imperfection. Finally, the coupling between 
the flexural modes plus the nonlinear terms in the equations 
of the axisymmetric mode were neglected. The resulting 


equations for the undamped motion were 


T 

e+ wo = Poe to (Sop 
mn 5 2 4 we 2 0 
wae pt y(n? ym C12" Cee n“we) oS ag 

7) 12 2 255 wi 
wl 1 2 T 
ee —) 
1 
f= pe eee (5 J6n 
Equation (5.6) can be written in the form 
n n 
W 2 Seen | Wee Za = 
Sr +n re w ) ial a ee ( 5 outed 
wl wi 
where 
2D ane 4 

pie ee ee (5.8) 

n 1 2 1 2 

(ne et =) 
l 
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The nondimensional static pressure which will cause classi- 
cal buckling of the ae mode is equal to Bee given» by (5.3). 
The natural frequency of the ae mode, Roe is obtained 
from (S97) by fequdaeing we to zero, from which 


0? = n“pm (5.9) 
nN Cr. 


Table II contains the natural frequencies, the natural 
periods Th = 27/2, and the buckling pressures of several 
of the modes of the shell under consideration. Note that 
the axisymmetric mode has the highest natural frequency, 
as wasS assumed in Chapter IV. Also, mode 6 has the low- 


est buckling pressure and, thus, is the static buckling 


mode .t 


Hyperbolic Response. Anderson and Lindberg [19] 


0 


pointed out that whenever w > ier the Se@ution of (5.77) 


is hyperbolic in nature, and the amplitude of the dis- 
placement grows in an exponential manner with time. The 


closed form solution of equation (5.5) is 


2 cc 


ee = 2 5 [e "To ~ cos t + = sin Tie (5.21805) 


1+T O 
Oo 





when A and ae are initially zero. As seen from Table II, 
the frequencies of modes four through eleven are consider- 


ably less than 2 the frequency of the uniform radial 


or 


mode. Consequently, these modes will not respond 


lMode 6 is the mode for which n = 6. In the sequel, 
the modes will be referred to in this manner. 


re 


TABLE II 


NATURAL FREQUENCIES, NATURAL PERIODS, 
AND BUCKLING PRESSURES PREDICTED BY 
EQUATIONS OF ANDERSON AND LINDBERG [19] 


n om Th pr x10" 
0 i C000 Ga 
Cees 2 45.48 lege 
5 0.1168 5 Saree eer 
6 0.1268 49.56 4.47 
7 Oza 6 1 40.25 Serre 
8 Q. 1 9)5ai oe. 20 5. OS 
G 0.2426 Zea face 
0 Ore 2 le 2 ieee So 
0 U. oe Lr. OW PO 
EZ 0.4260 14.85 12.43 
av Ur. 4g Zao 14.51 
14 0 .. Smee) Lor 9 6 L Gusts 
PS 0.6567 og. 36 Vo 7 


Sagal quantities are nondimensional. 


IFN 


appreciably to the harmonic oscillation of the axisymmetric 


motion. * Discarding the harmonic terms gives 


Pageea pala 


iT T 
Poto O 


w (tT) = ——e 


Since ees for the loadings of interest to the present 
study, the above equation reduces to 
0 


w (tT) = Poe 


which represents the quasi-static response of the axi- 
symmetric mode to the applied load. Substituting this 


result into (5.7) yields 


— + nate - pe ) — =n pe (Sri) 


As pointed out earlier, the solution of (5.11) is 


-T/T 
hyperbolic in nature when Poe espe Hence, aS long as 


the pressure ee exceeds the critical pressure of any 
one of the modes, that mode will grow in an exponential 
manner. Such modes will be referred to as "hyperbolic" 
modes. 

By virtue of their simplified analysis, Anderson and 
th 


Lindberg [19] were able to express the response of the n 


: coat : heer 
mode in terms of the amplification ratio w /wi , as can be 


-rhe Hesuleswor NMeoetvor and Lovell] wh2 Odueshew this to 
be true. 


IY, 


seen in (5.6). They showed that an amplification of 1000 
is tantamont to permanent buckling. The equations of 

the present analysis can not be solved directly for the 
amplification due to the nonlinear terms in the equations. 
However, once a solution is obtained, the amplification 
of each mode can be computed, based upon the values given 
to wi; and direct comparisons can be made with the 


results of Reference [19]. 


Mathieu Instability. A mode can also be excited to 
large amplitude through the mechanism of parametric 
resonance, or Mathieu instability. For Tort, the essen- 
tial features of the mechanism are present when (5.10) is 
reduced to 


-T/T 


oe = 328 (e ° - cos 1) (5 ie 


Substituting (5.12) into (5.7) and making use Of (5.9) sie 


EO 
n Z ae n 
w + (2 “n'P.e + Ti P, cos T) Ww = (5 See 
2 Tas n 
NP, (e = COs tT 
3 , 

it a and wi = Q, (5.13) reduces to Mathieu's equation 

wee (at + m1). cos tT) w = 0 (5.14) 


>The assumption that i 2a is not valid for the 
present problem. Consequently, the stability chart of 
(5.14) does not adeguately identify the Mathieu modes. 


nS 


There are two possible solutions to (5.14), both of them 
oscillatory. However, one of them grows unbounded with 
increasing time, while the other remains bounded. The 
unbounded growth is known as parametric resonance or 
dynamic instability. The name Mathieu mode, as used in 
the sequel, iS given to any mode that grows to large 
amplitude in an oscillatory manner. 

The stability boundaries for the solutions of (5.14) 


are well documented [7, 47] and are shown by Figure 5.1. 


2 
it a POs, n“p.) falls in the unstable region, the 


n 
solution for that mode will grow without bound. The 
Stability chart shows that when the frequency of a 

flexural mode iS approximately one-half of the axisymmetric 


natural frequency (2, = 1.0), the response of that flexural 


0 
mode will be unstable, even for a very small peak pressure. 
The points corresponding to modes 12, 13, 14, and 15 for 
the shell considered here are shown in Figure 5.1. 

Anderson and Lindberg [19] were interested in pre- 
dicting the critical buckling loads. From experimental 
results, they knew that dynamic buckling occurs as the 
result of very rapid growth of the buckling mode accom- 
panied by very little oscillatory motion. However, the 
Mathieu instability occurs in an oscillatory manner re- 
quiring several cycles before the displacements become 
large. Furthermore, when the amplitude of a Mathieu mode 


becomes large enough to cause plastic flow, energy will 


be dissipated, the axisymmetric oscillations will decay, 
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and subsequent growth of the Mathieu mode would most 
likely not occur.” Consequently, they were able to 


disregard the response of the Mathieu modes. 


McIvor and Lovell's Contribution 

Under impulsive loads, the initial growth of the 
flexural modes is governed by the Mathieu equation. The 
impulsive response of the axisymmetric mode can be ob- 
tained from (5.10) by taking the limit as TS goes to zero 


while the product Pala remains constant. The result is 
wo(T) =p T. sin T 153) 
o-G : 


Slibseteucing (5.15) mittee (5.7 jand making use of (5.9) 


lead to 


2 


n Z : Mae. lle? 
+ (82 -n pot, sin T)w = win Pot, sin t (orelie ) 


W 


If wi’ = 0, (5.16) reduces to the Mathieu equation 


n 2 


: 2 — 
w+ (2 meas Sin 1) weSa0 Csr? ) 


McIvor and Lovell's investigation [20] was concerned 
with the undamped response to purely impulsive loads. With 
an equation similar to (5.17), they identified the Mathieu 
modes from a stability chart. The total deflection was 
assumed to consist of the Mathieu modes plus the axi- 


symmetric mode. An energy balance revealed that the 


4 the axisymmetric oscillations are due to the initial 
conditions. Thus, there is a limited amount of energy in 
them. Once this energy is dissipated, they will no longer 
exist; and the mechanism for dynamic instability is removed. 


pez 


response of the neglected modes was, in fact, negligible. 
All of the nonlinear coupling between the modes was re- 
tained. During the response, the modes underwent beat-like 
oscillations accompanied by a cyclic exchange of energy 
between the modes. The maximum amplitude of a mode did 

not always occur during the first beat, sometimes occurring 
during the second beat. The stresses developed during the 
response to an impulse were several times higher than the 
stresses of an axisymmetric response to the same impulse. 
Furthermore, the largest stress usually did not occur 
during the first Significant exchange of energy, but 
occurred later when several modes were nearly in phase. 

The maximum amplitudes attained by the flexural modes were 
nearly independent of the magnitudes of the perturbation 

in each mode; however, the phasing between the modes was 


affected by the magnitudes of the perturbation. 


Ifil. INVESTIGATION OFSUEE PROBE EMS lime 


METHODS. DEV ERO ER De Niet Hilo oma Opn 


Preliminaries 


Minimum Number of Axial Grid Points. In order to 
minimize the computing time, the number of axial stations 
was kept as small as reasonable accuracy permitted. The 
criteria used for assessing the adequacy of a given number 
of stations were: (1) the accuracy of the circumferential 


in-plane forces, (2) the accuracy of the amplification 


ee 


factors, and (3) the accuracy of the representation of the 
axial variation of the radial displacements. 

Since the load and the initial imperfection were 
assumed to be symmetric with respect to the midspan of the 
shell, it would not be unreasonable to expect the response 
to also be symmetric with respect to the midspan. If true, 
this would effect a considerable savings in the computing 
time. Two fully nonlinear solutions were computed to 
determine whether or not the response was indeed symmetric. 
The first run employed the simple support conditons at 
both ends of the shell and fourteen axial stations. The 
second run used the midspan symmetry conditons and seven 
axial stations. The displacements of the two solutions 
were identical to at least four significant figures. 
Therefore, the response for this problem was assumed to be 
symmetric with respect to the midspan of the shell. 

A series of computer runs was made using a different 
number of axial stations over the half-length of the shell 
for each run. There were no significant differences in 
the results when nine, eleven, or thirteen axial stations 
were used along one-half of the shell. However, when the 
number of stations was reduced to seven, the error in the 
circumferential in-plane force was approximately 7 per cent 
near the end of shell; and smaller errors were noted in 


the radial displacements and the amplification factors. 


2 3 


Consequently, nine stations over the half-span of the shell 
were used for the solutions presented in the following 


sections of this chapter. 


Selection of the Participating Modes. Since the 
computer program is dimensioned to accomodate a maximum 
of fifteen Fourier od ooe the series expansions for the 
displacements had to be truncated. However, the truncated 
series was not restricted to the first fifteen modes, 
modes 0 through 14. Instead, the participating modes 
were selected to include only those expected to respond 
to the particular loading condition. The procedure used 
for selecting the responding modes can best be explained 
by an example. 


Consider the loading 


-t/ 2008 


P(t) = 70 e (S alae 


Comparing the 70 psi (Dy =O ioe peak pressure with 


the buckling pressures listed in Table II reveals that modes 
5 through 8 should exhibit hyperbolic behavior during the 
early phases of the response. Hence, these modes were 
included in the series. 

Identifying the Mathieu modes is a little more diffi- 
cult. As a result of the assumptions made to arrive at 


(5.14), the stability chart of Figure 5.1 is not Suftitilem@ege 


>The number of modes is limited primarily by com- 
puting time. 
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to identify the Mathieu modes. Conseguently, a procedure 
was used that investigated the stability of each of the 
possible Mathieu modes. For example, Table II shows that 
modes 12 through 15 are possible Mathieu modes since they 
have frequencies iying in the neighborhood of one-half 

the frequency of the axisymmetric mode. To determine 
which of these modes were Mathieu modes, the series was 
truncated so that only modes 0 and 12 through 15 were 
retained. As a result of this choice of Fourier indices, 
each of the flexural modes couples only with mode zero.° 
Finally, the coupling terms in the equations of mode 0 
were equated to zero by a simple alteration of the computer 
program. This eliminated the mechanism for the extraction 
of energy from the mode zero oscillation. The resulting 
eguations of motion of the flexural modes resembled 
Mathieu's equation, and any mode that grew to a large 
amplitude under these circumstances was assumed to be a 
Mathieu mode. A computer run was carried out based on 


this procedure using the pressure given by (5.18). The 


©the modes which couple are determined by forming all 
combinations of sums and differences of the Fourier numbers 
in the series. If any sum or difference equals any of the 
Fourier numbers in the series, then the two Fourier numbers 
contributing to that sum or difference are said to couple 
with the third mode. For example, modes 3 and 5 couple to 
influence modes 2 and 8 since 3 + 5 = 8 and 5 - 3 = 2. Mode 
13 couples with itself to influence modes 0 and 26 since 
13 - 13 = 0 and 13 + 13 = 26. For a more detailed explana- 
tion see Appendix B, Reference [30]. 
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results are given in Table III. The maximum amplification 
listed in the table is defined as the maximum value of the 


ratio 


=— (5. ie 


where station i is the whole station nearest to the midspan 
of the shell. From Table III, mode 14 can be readily 
identified as a Mathieu mode. However, as shown in Figure 
5.1, the Mathieu stability chart of (5.14) indicates scm 
mode 13 should be the unstable mode. 

A possible explanation of this discrepancy can be 
surmised by considering (5.13). When Ws equal to zero, 
(5.13) reduces to 


-t/T 
‘ae (ae - n“p.e ai np. cos 7a = 0 (52100 


TABLE LIL 
RESPONSE OF INDIVIDUAL MODES WITHOUT 


EXTRACTION OF ENERGY FROM AXISYMMETRIC MODE 


Dic 6007 —x MWnes eee 20k 
O O 


Maximum Amplification 
Mode Attained During Interval 
tT = 0 to 140 


ay eae 
cig Ges 
14 0). 
15 3.4 
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This equation is Similar to the Mathieu equation, except 


that 0% has been replaced by 


Be a te (5.21) 
= oe 
nn © 


The quantity expressed by (5.21) can be thought of as a 
time varying frequency of the re mede. if this quantity 
is plotted along the abscissa of the Mathieu stability 
@m@art instead of cee as in Figure 5.2, the points associated 
with the flexural modes 12, 13, 14 and 15 move horizontally 
across the chart as Tt increases. 

Examination of Figure 5.2 reveals that at Tt = 0, 
modes 12 and 13 are in the stable regions, and modes 14 and 
15 are unstable. In computing the amplifications listed 
in Table III, the integration was terminated at t = 140. 
During that time interval, modes 12 and 13 grew very little; 
mode 14 was amplified nearly two thousand times; and mode 
15 was amplified 3.4 times. These results can be correlated 
with the trajectories of the corresponding points in Figure 
5.2. During the interval of integration, modes 12 and 13 
remained in a stable region, but mode 14 remained in the 
unstable region. Mode 15, which was initially in an un- 
stable region, quickly moved into a stable region. Since 
growth to a large amplitude requires several oscillations 
of the axisymmetric mode, mode 15 moved into a stable region 
before having time to grow appreciably. Had the integration 
continued, mode 13 would have entered the unstable region; 


but mode 12 would have remained stable. 
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MATHIEU STABILITY CHART FOR EXPONENTIALLY 


DECAYING UNIFORM PRESSURE 


; ae 4, | , 
PB =70 psil(p.=6.3 7x10 "), t=3.0msec (t= 208) 
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As a consequence of this correlation the stability 
chart of Figure 5.2, which is based on the simplified 
equations of Reference [19], was considered to be adequate 
for identifying the Mathieu modes for the problem under 
consideration. To identify the Mathieu modes for a prob- 
lem with different loading conditions, it might be neces- 
sary to make several preliminary computer runs similar to 
the one described above. 

When computing a fully coupled nonlinear solution 
mor the loading given by (5.18), modes 0, 1, 4, 5, 6, 7, 
8, 9, 10, 13, 14, and 15 were retained in the Fourier 
expansions. Note that modes 4, 9, and 10, modes in the 
neighborhood of the hyperbolic modes, were retained. The 
Mathieu modes that were retained included all modes whose 
trajectories in Figure 5.2 passed through the unstable 
region. Finally, mode 1 was also included because it be- 
comes involved in the coupling of all the modes whose 
wave numbers differ by sake ate, The adequacy of this choice 
of modes was established by computing a solution in which 
modes 0 through 14 were retained. The extra modes, 2, 3, 
ll, and 12, did not respond appreciably; and including them 


did not alter the response of the other modes. 


Tin the Donnell theory mode 1 is inaccurately 
described for many classes of shells. However, if the shell 
is thin and L/as2, as is the case here, the representation 
of mode 1 is reasonably accurate [48]. 
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Results 
Computer solutions were generated for four specific 
cases, hereafter referred to as A, B, C, and a In cases 


A and B, the damping coefficient c was taken equal to zero, 


while in cases C and D, c = 3.416 x ow >: Cases A and C 
were each composed of three runs with ty = 3.0 msec 

(TO = 208) and Po equal to 60 psi (Do = 5. 460 107°), 70 
psi (DO =i6.. 37 es iowa and 77 psi (Do =a Olle ikon 


Cases B and D were also composed of three runs each, but 
at three different peak pressures with to = 1.0 msec 

(TO = 69.3). The results for these cases will be compared 
with the results obtained by Anderson and Lindberg [19] 
and McIvor and Lovell [20]. 

The significance of the amplitudes of the initial 
imperfection modes was also investigated by computing a 
solution in which the amplitudes wi were taken from a 
table of random numbers. Finally, some results are pre- 
sented to show the importance of the nonlinear coupling 


between the modes 


Case A. The results for this serves et. tna 
summarized in Table IV. The first column shows the Foun 
indices of the modes included in the solution. The seco 


column shows the maximum amplification of each flexural 


mode, as given by (5.19). The third column contains the 


Ba typical solution required approximately 45 minutes 
on the IBM 360/67. 


0) 


TABLE IV 


CASE A, @MPLIFICATIONS AND STRESSES 


— = 3.0 msec Tas = 208), ZERO DAMPING 
Bue psi : = Shtdex wer 
O P 7 Po , 
Maximum circumferential stress = 21,600 psi 
Nexium circumferential stress occurred at T = 116 
Integration terminated at T = 118 
mode maximum amp. t fOr Maxam aie. 
0 
|| 72 aa0.4 
5 a? 2 52 .6¢ 
6 AO TL.05 
7 PARAS 52.78 
8 9.18 Se 
9 4.60 Dey, 13 
ies 1,94 LihLg.6 
14 97.39 De. 
es IL Aes 16.6 
P= 70 psi - = 6s eeu 
aS p 2 Ee ° 
Maximum circumferential stress = 32,200 psi 
Mem imum circumferential stress occurred at +t = 165 
Integration terminated at Tt = 198 
mode maximum amp. tT for maximum amp. 
0 
il FO 8 190.6 
4 pee. 88.35 
5 26. 7@ yu. 7k 
6 222. 3 97.35 
7 ae 7 We .23 
8 22 Oe AD. 37 
9 7.84 26.41 
"6 4 le5 19.43 
is ano 3.8 
14 Zoe? 164.6 
15 4.06 162.8 


id. 


TABLE Lv  “COneLlmucs 


= -4 
ie = 77 psl ee Te Onl Ee) 
Maximum circumferential stress = 64,300 psi 
Maximum Clreumierent lal estress™oCcecuLuCa cee t]—=— . 
Integration terminated at t = 166 

mode maximum amp. t for maximum amp. 
0 
ic L268 Lo a 
4 LO. V5 Sre7 
5 esta 7 3 oNeO 
6 Laan P24 
7 2586 90.96 
8 97220 Mis 
9 25 oO 118.9 
10 7.79 lease 3 
is 116.6 ira So 
14 2395.5 1Aay i 
LS Ae a ewe 


nondimensional time at which the maximum amplification 
occurred. The magnitude of the maximum circumferential 
stress at @ = je as given by (4.31), and the nondimen- 
Sional time of its occurrence are also listed. Finally, 
the nondimensional time at which tne run was terminated 
is indicated.” 

The time histories of the radial displacements at 
the whole station nearest to the midspan are shown in 
Figures 5.3-5.5 for the four or five modes that attained 


the largest amplifications. In addition, the Figures are 


? the nondimensional period of the axisymmetric mode 
1s approximately equal to 27. 
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annotated with values of the circumferential stress to 
show the times at which a relative maximum occurred. Note 
the different scales for w’ in the three Figures. 

From Table IV, it can be seen that when the peak 
pressure was 60 ps1 (Do = 520 Cr ae hyperbolic modes 
5, 6, and 7 grew to large amplitudes, and that the growth 
occurred during the early phases of the response. Exam- 
ination of Table II shows that these are precisely the 
modes having static buckling pressures less than the peak 
pressure of the load pulse. Furthermore, mode 6, the 
static buckling mode, attained a larger amplification than 
the other hyperbolic modes. If comparisons are made for 


the other two runs, similar observations can be made. 


In the third run, mode 1 grew to an amplification of 
126.8 just prior to the peaking out of the hyperbolic modes. 
This growth of mode 1 was very likely due to coupling 
between the hyperbolic modes that produced large coupling 


terms in the mode 1 equations. The reason for including 
mode 1] @nwethemtruncated sericemrs now clear. 


In all three runs, mode 14 attained the largest am- 
plification of the Mathieu modes. Inspection of Figures 
5.3-5.5 reveals that mode 14 underwent beat-like oscilla- 


tions requiring several cycles to attain a maximum ampli- 


tude. However, it was not clear in any run that the largest 


amplitude of the Mathieu modes had occurred prior to the 


termination of the computer run. In fact, the information 
shown in Figure 5.2 indicates that mode 13 would have grown 


large for an. 
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Comparison of the time histories for the three runs 
reveals that at the lower peak pressures, the Mathieu 
mode dominated the response. As shown in Figure 5.3, the 
amplification attained by mode 14 was more than twice as 
large as the amplification of any of the hyperbolic modes. 
Furthermore, the largest stresses occurred as a result of 
the growth of the Mathieu mode. As the peak pressure and 
the total impulse Pot, were increased, the amplifications 
of both the Mathieu mode and hyperbolic modes were of 
about equal magnitude, as shown in Figure 5.4. Figure 
5.5 shows that at the highest pressure and impulse, the 
hyperbolic modes dominated, and the maximum stresses 
occurred as a result of their growth. 

In Figure 5.6, the maximum amplification of the hyper- 
bolic mode 6 is plotted against the peak pressure and, 
equivalently, the total impulse. Assuming an amplification 
of 1000 to be the dynamic buckling ih ea Rugiare ~5...6 
shows that the buckling pressure is 76 psi when the total 
impulse is 228 psi-msec. Note the logarithmic scale for 


the amplification factor. 


Case B. For the runs of case B, the time constant 
of the decaying load was reduced to 1.0 msec (T = 69.3); 


and the pressures for the three runs were 100 psi 


(p, = 9.10 x 10°"), 110 psi (p, = 10.01 x 10°"), and 120 


psi (Do OR oS 2 o5e we The hyperbolic modes were 5, 6, 


LOS Dani and Lindberg [19] also used this criterion. 
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FIGURE 5.6 
CASE A 
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MODES 
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7, 8, and 9, all of which have buckling pressures less 
than the peak pressures of the loads. In addition, modes 
4 and 10 were included. The Mathieu modes were iden- 
Pitted freMmea Stability Chart similar to Figure 5.2 and 
included modes 13, 14, and 15 at the lowest peak pressure 


plus mode 12 at the two higher peak pressures. As with 


case A, the damping coefficient was set to zero. Two of 
the runs were computed out to * = 200; and the third, to 
T = 150. The results are summarized in Table V. 


As in Case A, the character of the response changed 
considerably as the peak pressure and the total impulse 
were increased. When Py = 100 psi (Do = 9.10 x Ome) 
the maximum amplifications of the Mathieu modes were 
roughly twice as large as the amplifications of the 
hyperbolic modes. For the second solution, the peak pres- 
sure was equal to 110 psi; and the maximum amplifications 
of the hyperbolic and Mathieu modes were approximately 
equal. The loading of the final run was near the critical 
dynamic buckling limit and boosted the hyperbolic mode 6 
eewan amplificatzon of 832.6. The largest amplification 
of the Mathieu modes was 340.6. The maximum stress occurred 
well after the hyperbolic modes had passed their peak 
Bioeth cations @n Gre first and second reruns. On the other 


hand, the largest stress occurred in conjunction with the 


large response of the hyperbolic modes in the last run. 


ioe) 


TABLE V 


CASE 3B, AMPLIDICAT TONS. 2) (pe ot Bea. 


ty = 1.0 msec as = 69.3), ZERO DAMPING 
p= Gast = oe 10 
3 p ae ° 
Maximum circumferential stress = 55,000 psi 
Maximum circumferential stress occurred at T = 199 
Integration terminated at T = 200 
mode maximum amp. t for maximum amp. 
0 
1 ere ss 1S yaar 
4 S gined) RSS) «6 
5 Bron 5 Ore 
6 ee? 65.44 
i Loy Bee 
8 Sula, 40.24 
2 A038 ilgsyera | 
a) Oe Oar 
ie 3 os 182.9 
ies Jee co Gene 
14 2a 550 7 Oe 
is Jeo ed ORO 
- ; . _~ 4 - -4 
Ee = 110 psi = Py JOO sia. eG 
Maximum circumferential stress = 65,500 psi 
Maximum circumferential stress occurred at T = 117 
Integration terminated at T = 200 
mode maximum amp. t for maximum amp. 
0 
i 0 oe ou 
4 14.82 L7 OG 
5 Be. 9 Gis. 516 
6 291.5 Jie 
i 2 Oem Smee 
8 ore. 1 46.52 
5) 54.66 ees 
10 ie cic 3 led 
ie Slee. 6 1956 
12 22 oes 184.2 
Ne 2 ioe 184.1 
14 ZS 2 5 
His) L286 Eo ae 
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TABLE vo eecentinued 


_ -4 
es = 120 psi : Po = G92 x LO 
Maximum circumferential stress = 76,200 pee 
Maximum circumferential stress occurred at tT = 57 
Ingegration terminated ater = 150 

mode maximum amp. t for maximum amp. 
0 
1 148.6 oe SL 
3 36.77 OO... 
4 34.58 Le. 
5 eG. 5 Ble 4 
6 832.6 Bio) 
y 656.4 66.93 
8 218.2 527006 
9 FO. 54 7.54 
10 25s 28204 
ek 30.94 35°. 7 
12 152.8 ied 
13 pea a ee 
14 34 0.6 ielesi6 
cS 250.3 98.83 





The behavior of the Mathieu modes reflected the 
shorter time constant and the higher peak pressures. In 
all three runs, all of the included Mathieu modes grew to 
Significant amplitudes. Again this can be explained by 
emo cadoility chant “of "Fraquremer 2 .c.As n“p. 1s increased, 
the trajectories of the Mathieu modes are shifted upward 
and extended in length to the left; the reduced time con- 
SUdiltemeauses the™points corresponding to the modes to 
move across the chart at a faster rate. The net result 
Mdcmeladt MOdeS I2,0o dnce | Seuas Well as 14) were: immitehne 


unstable region sufficiently long to respond. 
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Note that in all three runs, the stress exceeded 42 
ksi, the yield point stress of the shell material, while 
the amplifications of the hyperbolic modes never exceeded 
1000, the dynamic buckling limit used by Anderson and Lind- 
berg [19]. 

The maximum amplification of the hyperbolic modes as 
a function of the peak pressure and the total impulse is 
shown in Figure 5.7. An amplification of 1000 corresponds 


to a pressure of approximately 122 psi. 


Comparison with Anderson and Lindberg's Findings. 


Anderson and Lindberg [19] identified the peak pressure 

P and the total impulse Pot, as the significant parameters 
in the dynamic buckling problem. Their dynamic buckling 
curve, along which the maximum amplification of the hyper- 
bolic modes was 1000, is shown in Figure 5.8. The sharp 
break in the curve divides the regime of the elastic model, 
on the lower right, from the regime of the tangent modulus 
model, at the upper left. As the impulse of theload goes 
to infinity, corresponding to a step load, the dynamic 
buckling curve approaches the static buckling pressure 
line, The dynamic curve shows that the shell can endure 

a pulse having a peak pressure larger than the- static 
buckling pressure, with the critical peak pressure depend- 
ing on the total impulse of the load. While Anderson and 
Lindberg's model gives the dynamic buckling curve, it does 


not predict the response in the region lying below the 
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dynamic buckling curve since the Mathieu modes were not 
included. However, the algorithm of this dissertation is 
capable of computing the displacements and the stresses 
in the subcritical region, notably, in the region between 
the static buckling line and the dynamic buckling curve. 
The two straight lines running diagonally across 
the Figure are lines along which ce 1s equal to 3.0 msec, 
case A, and 1.0 msec, case B. The two points labeled 
CASE A and CASE B are the buckling pressures and impulses 
which were obtained from Figures 5.6 and 5.7. They are 
based on an amplification of 1000. Note that the points 
fall very close to Anderson and Lindberg's buckling cruve. 
Therefore, Anderson and Lindberg's treatment of the boundary 
conditions and their neglecting of certain nonlinear 
coupling terms do not lead to any Signigicant errors in 


Predicting the amplification of 1000. 


Comparison with McIvor and Lovell's Findings. Under 
an impulsive load, as considered by McIvor and Lovell [20], 
the axisymmetric mode oscillates about its unloaded static 
equilibrium state, and the quasi-static response of the 
axisymmetric mode is zero. Therefore, the hyperbolic 
modes, which respond primarily to the quasi-static motion 
of the axisymmetric, are not excited. In cases A and B, 
the amplitudes of the hyperbolic modes were small at the 
lower pressures and impulses, because the quasi-static 


response of the asSisymmetric mode was very small; yet the 


145 


high frequency oscillation of the axisymmetric mode was 
sufficient to excite the Mathieu modes. The observed 
behavior of the Mathieu modes agrees qualitively with 
McIvor and Lovell's findings.** For instance, in Figure 
5.4, the Mathieu mode responded in a beat-like manner 
requiring several cycles to attain the first peak. Fur- 
thermore, the peak amplitude of the Mathieu mode was 
larger in the second beat than in the first, leading to a 


large stress. 


Maximum Stresses. Although the theory of Reference 
[19] appears to be very adequate for computing buckling 
thresholds, the results for cases A and B, contained in 
Tables IV and V, indicate that the stresses may exceed 
the yield point stress at peak pressures and total impulses 
well below the buckling threshold. When the load was sub- 
critical, maximum stress occurred well after the lower 
frequency, buckling modes had attained their maximum 
amplitudes. However, at loads near critical the maximum 
stress occurred at about the time the dominant hyperbolic 
mode was at its largest amplitude. 

Since a structure designed for repeated use, such 
as the NWL shock tube, would experience fatigue weakening 
if the stresses exceed the yield point stress, one might 


ask whether a normal amount of viscous damping would 


aa quantitative comparison can not be made due to 
the differences in the analyses. 
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attenuate the Mathieu response, thus reducing the stresses. 


This aspect is considered in cases C and D. 


Case C. The loading of case C was the same as the 
loading of case A except that viscous damping was included. 
The nondimensional coefficient c was taken equal to 
Be416 x Ome The maximum amplification of each mode and 
the time of itS occurrence, the maximum circumferential 
stress and the time of its occurrence, and the time at which 
the integration was terminated are contained in Table VI. 
Pieeaddition, Figures 5:9-5./ienow the time histories of 
the active modes, annotated to indicate the times at which 
large stresses occurred. 

The computer runs of case C were carried out for 
longer periods of time than were the runs of case A in 
order to include any Significant growth of mode 13 that 
might lead to high stresses. The amplifications of 
Table VI indicate that mode 13 had moved into the unstable 
region of the Mathieu chart of Figure 5.2 prior to the end 
of each run. On the first two runs, mode 13 became un- 
Stable only after the damping effects had extracted much 
of the energy from the motion of the shell. Consequently, 
mode 13 did not grow sufficiently to generate a higher 
Stress than had previously occurred. Due to the higher 
peak pressure of the third run, mode 13 became unstable 


much sooner. However, the mode 13 response was not respon- 


Sible for the maximum stresSs Since the motion of the 
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TABI 7 


CASE C, AMPLIFICATIONS AND STRESSES 


ce = 3.0 msec (Us = 208), DAMPING INCLUDED 
ERGO. ped ee eae 
O P 7 By ° 
Maximum circumferential stress = 16,250 
Maximum circumferential stress occurred at T = 136 
Integration terminated at T = 300 
mode maximum amp. tT for maximum amp. 
0 
il 94 287.5 
4 1.49 21 65 
5 iOeiG k Dos 
6 38.24 20. 
q 22.98 52.73 
8 Seo Ss 33.86 
9 Ae 202 
13 22.86 296.0 
alee Tiere ey 
15 ie 45 84.41 
2 = 7G joc : Seevaeeo 
O P 7 Py : 
Maximum circumferential stress = 26,400 psi 
Maximum circumferential stress occurred at Tt = 110 
Integration terminated at Tt = 257 
mode maximum amp. Tt £Or maximum dip. 
0 
is A 56 S5. 72 
4 1.99 3 3 ow 
Ss 24279 71.78 
6 202720 97.66 
7 OSes 72.14 
8 Bilee2iG 46.44 
9 7.68 26es4 
16 4.09 19.41 
13 4.98 252 
ina Aa ep? a3 
IL U9 Jem 32 
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TABLE VI Continued 


a ; 2 -4 
Py = 77 psl eee loa x 0 
Maximum circumferential stress = 53,800 psi 
Maximum circumferential stress occurred at Tt = 110 
Integration terminated at tT = 198 
mode maximum amp. t for maximum amp. 
0 
1 91.06 iOS 9 
4 3s IESG 5 
5 87.00 1976 
6 1S At her ae 
q 396.2 91.54 
8 69.72 eee 
9 eye alt 119.0 
10 6.54 168.1 
13 IS Se 30 197.0 
14 di 4 LO325 
ib5) 9.92 3 30 


Mathieu modes was greatly overshadowed by the early growth 
of the hyperbolic modes. 

The effectiveness of the damping in attenuating the 
Mathieu modes can be ascertained by comparing the mode 14 
amplifications of the damped and undamped responses. 
Comparison of the amplifications of mode 14 contained in 
Table IV with those contained in Table VI reveals that the 
maximum amplification of mode 14 was reduced by 25 to 30 
per cent by the inclusion of damping. 

Assuming that the damping would have extracted 
sufficient energy to preclude the development of stresses 
higher than those computed in Table VI, a critical pressure 


and impulse based on the yield point stress can be 


a 2 


determined. In Figure 5.12 the maximum circumferential 
stress is shown as a function of the peak pressure. The 
critical pressure and impulse corresponding to a maximum 
stress of 42 ksi are 75 psi and 225 psi-msec, respectively. 
Note that the yield stress occurs at a pressure where the 
hyperbolic modes have much larger amplifications than the 
Mathieu modes. However, the bending stress, given by the 
second term of (4.29), depends on aC. consequently, the 
Mathieu modes with their higher mode numbers can make 
Significant contributions to the stress even though their 


amplitudes may be small. 


Case D. In case D, the time constant Ps was taken 
equal to 1.0 msec (ity = 69.3), the time constant used in 
case B. The nondimensional damping coefficient was again 
taken equal to 3.416 x iG iier and three computer runs were 
carried out at peak pressures of 90 psi, 100 psi, and 110 
psi. The results are summarized in Table VII. Note that 
for all three runs, the integration was carried out long 
enough to include the maximum amplification of mode thir- 
Teen. 

Comparison of the amplifications of the first run of 
case B (Table V) with the second run of case D (Table VII) 
shows that the amplifications of the Mathieu modes were 
again reduced approximately 25 to 35 per cent. However, a 
Similar comparison of the third run results of case D with 


the second run results of case B reveals that at the higher 
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FIGURE 5.12 
CASE C 


MAXIMUM CIRCUMFERENTIAL STRESS vs. 
PEAK PRESSURE 
tp = 3.0 msec (T = 208), Damping Included 
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TABLE VII 


CASE D, AMPLIFICATIONS AND STRESSES 


ee = 1.0 msec oo = 69.3), DAMPING INCLUDED 
_ z - -4 
es = 90 ps1 > Dy Spero x 1.0 
Maximum circumferential stress = 32,100 psi 
Maximum circumferential stress occurred at T = 187 
Integration terminated at T = 224 
mode maximum amp. tT for maximum amp. 
0 
i. 26 ak 2: DAS 
4 2.44 237 
5 19.15 5259 
6 Gio 4 Boieo 5 
7 53a 46.70 
8 24.94 40.03 
9 Tie 76 26.39 
10 6.47 19.56 
ik 4.66 DAG 
es 199.9 5322 
14 15 5405 89.83 
5 OS 49.03 
> SLO) joe oe 90x anes 
Oo 246) : 
Maximum circumferential stress = 37,900 psi 
Maximum circumferential stress occurred at T = 182 
Integration terminated at T = 193 
mode maximum amp. T for maximum amp. 
0 
1 S767 yee 2 
4 374 G2. 2 
5 20% 2 58.60 
6 12506 65. 44 
7 ies. 6 S270 6 
8 50.208. 40.28 
9 18.98 26.93 
10 9.46 19.84 
ek Gai rez ea 
3 2ioieS Me? . 5 
14 208.3 83.39 
15 S426 7 56.14 
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TABLE VII Continued 


= . = -4 
PS = 110 psi aos OO isc 80 
Maximum circumferential stress = 45,400 
Maximum circumferential stress occurred at tT = 120 
Integration terminated at tT = 188 
mode maximum amp. tT for maximum amp. 
@) 
il 56.98 4 AO 
4 Toe 1359 
5 52.84 Goma i 
6 267.9 Alco 
d 2450 59.19 
8 LOBral 46.47 
9 3300 sei 
10 oa A0R0S 
Pa 9.88 E7655 
eZ 68.58 183.6 
ile 2 Sujal 146.5 
14 252 P1828 
RS" EO Os 5 62.552 


peak pressure damping was not nearly so effective in re- 
ducing the amplifications of the Mathieu modes. 

In Figure 5.13 the maximum circumferential stress is 
plotted against the peak pressure. Assuming that the yield 
point stress is 42 ksi, the critical peak pressure is 
approximately 106.5 psi. At this pressure the amplitudes 
of both types of modes are approximately the same. Thus 
the Mathieu modes play a dominant part in the stress state 


of the cylindergror this loading @eendiction. 


Critical Loading Based on Maximum Stress. A critical 


loading curve based on the maximum circumferential stress 
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FIGURE 5.13 
CASE D 


MAXIMUM CIRCUMFERENTIAL STRESS _ vs. 
t=1.0 msec ie = 69.3), Damping Included 


at 68 = 0° being equal to the yield point stress is plotted 
in the pressure-impulse plane in Figure 5.14. The curve 
for an amplification of 1000 obtained by Anderson and 
Lindberg [19] is also shown. Although cases C and D 
provide only two points to define the curve, clearly, 
critical stresses occur at pressures and impulses less than 


those required to produce amplifications of 1000. 


Axial Wave Form of Responding Modes. Anderson and 
Lindberg [19] assumed the axial variation of the radial 
displacements of the flexural modes and the initial imper- 
fection modes to be a half sine wave. The axial variation 
of the displacements was completely unspecified in the 
present analysis. For the solutions of this chapter, the 
initial imperfection modes were assumed to vary as a half 
Sine wave in the axial direction, corresponding with 
Anderson and Lindberg. In Figure 5.15, the axial distri- 
butions of the computed radial displacements in modes 6 
and 14 are plotted at arbitrarily chosen times. Note that 
the axial variation of mode 6 is essentially a sine wave 
as assumed by Anderson and Lindberg. On the other hand, 
the radial displacement of the Mathieu mode is not a 
Single half sine wave. 

The axial variation of the Mathieu mode shown in 
Figure 5.15 appears to contain a component of the form 
sin(37&/1) as well as a component of the form sin(7Té/1). 
The natural frequency 0 of a mode of the form 
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© n=6, T= Jel 
— half sine wave 
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wan sin (S25) cos nég 


is given by 
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A few quick calculations show that when gq is equal to 


three, the natural frequencies are equal to 


Te 23n 
10 Hee ei 
ie) ~439 i) 
12 ~-4964 
13 -5630 
14 moo] S 
The frequency for n = 12 is very close to one-half, and 


the possibility exists that this mode might be para- 
metrically excited by the oscillation of the axisymmetric 


mode, 74 


Rondel lmpewicecie sens. MEnmDEaActICce, the Claaneer 
does not know the magnitude of the initial imperfection 
until after the shell has been constructed. This leads 
one to examine the sensitivity of the solution to the magni- 
tudes of the initial imperfections. One computer run was 


made at P ee On Silies = 3 msec, and wi? Sol 0x Ouse 


teMcIvor and Lovell [20] recognized this possibility 


and included similar Mathieu modes in their calculations. 
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The magnitude of the imperfection in each mode was selected 
from a table of random numbers. The results obtained for 
uniform imperfections and random imperfections are shown in 
Table VIII for purposes of comparison. 

Note that the amplifications of modes four through 
ten are nearly the same in both cases. However, the dis- 
placements are Significantly different. On the other hand, 
the magnitude of the displacement, instead of the amplifi- 
cation, remainS approximately constant for modes 13, 14, 
and 15. Thus, it appears that the hyperbolic modes are 
very sensitive to the magnitude of the initial imperfection, 
whereas the response of the Mathieu modes is largely in- 
dependent of the initial imperfection. This is in agree- 
ment with the conclusions of Anderson and Lindberg [19] 
and McIvor and Lovell [20]. Anderson and Lindberg con- 
Sidered the hyperbolic modes to be directly proportional 
to the imperfection; while McIvor and Lovell found that 
the maximum amplitude of a Mathieu mode was rather in- 
sensitive to the initial imperfection. 

Since the stresses depend on the magnitude of the 
displacements, the stresses are very sensitive to the 
magnitude of the initial imperfection when the hyperbolic 


modes play a dominant role. 


Significance of Intermodal Coupling of the Flexural 
Modes. The maximum amplifications of the hyperbolic modes 


obtained by Anderson and Lindberg [19] agree well with the 
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amplifications computed by the fully coupled nonlinear 
theory. Since their analysis neglects, among other things, 
the intermodal coupling of the flexural modes, the possi- 
bility arises that the coupling between the flexural 

modes is of negligible importance. Neglecting intermodal 
coupling of the flexural modes is equivalent to retaining 
only those nonlinear products which contain the Fourier 
coefficient of the axisymmetric mode. Motion in the 
axisymmetric mode affects the motion of each of the other 
modes, and conversely each mode couples back to the axi- 
symmetric mode. However, modes i and j do not affect each 
other if neither i nor j is equal to zero. 

The computer program can be altered very readily to 
effect this limited coupling. A run was made for the 
parameters of case A, Po = 70 psi, using the altered 
program. The maximum amplifications for both complete 
and limited nonlinear coupling are presented in Table IX. 
The integration was terminated at t = 200 so that mode 
13 did not have time to grow. The magnitudes and times 
of the maximum amplification of all of the modes are nearly 
the same for both cases. On the basis of this one solu- 
tion, it appears that the intermodal coupling of the 
flexural modes for this problem is quite weak and can be 
neglected. 

Additional computer runs should be made in the future 
to cover a larger spectrum of loading conditons. If it can 


be shown that the intermodal coupling of the flexural modes 
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TABLE IX 


AMPIIE LT CATIONS WITH TOTAL 


AND LIMITED INTERMODAL COUPLING 


ae = 70 psi, Ge = 3.0 msec, wi = LOX ioe 


ZERO DAMPING 


TG eee @ elaine Limited Coupling 
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Maximum T FO Maximum tO 
Mode | Amplification | Max. Amp.| Amplification | Max. Amp. 


0 
1 7.933 | 190.6 039 198.1 
4 2.220 88.35 1.974 33.53 
5 26.70 palegia 26.20 71.78 
6 222.3 )) Gees 219.6 97.48 
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is indeed negligible for all loadings of interest, the 
analysis and the computer program can be simplified very 


Substantially. 


Significance of Coupling Terms in Mode Zero Equations. 
The coupling terms in the equations of mode 0, the axi- 
symmetric mode, provide the mechanism for the extraction 
of energy from that mode. A natural question to ask is 
whether the response of the flexural modes would be 
affected if this mechanism were removed. If this coupling, 
aS well as the intermodal coupling of the flexural modes, 
is unimportant, then the response of all of the flexural 
modes could be computed one mode at a time, using equations 
(5.5) and (5.6), or similar equations. Three computer 
runs were made to study this guestion. In all three runs, 
the coupling terms in the mode 0 equations were equated 
to zero. Run one included the flexural modes 5 through 8, 
run two included the flexural modes 12 through 15, and run 
three included the flexural modes 16 through 18. The 
choice of these flexural modes rules out any intermodal 
coupling of the flexural modes. The loading had a peak 
pressure of 70 psi and a time constant of 3.0 msec, the 
same as the loading of the second run of case A. The 
Maximum amplification of each mode from both solutions is 
contained in Table X. Comparison of the appropriate 
amplifications from Table IV with those contained in Table 


X reveals that the amplifications of the hyperbolic modes, 
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TABLE X 


RESPONSE OF INDIVIDUAL MODES WITHOUT 


ENERGY EXTRACTION FROM THE AXISYMMETRIC MODE 
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modes 5 through 8, were not changed appreciably when the 
coupling terms in the mode 0 equation were neglected. On 
the other hand, the amplification of the Mathieu mode, 
mode 14, was drastically larger. This is reasonable since 
the initial conditons rather than the loading cause the 
oscillations of the axisymmetric mode. Thus, the oscilla- 
tory energy of the axisymmetric mode is limited. If this 
limited energy can not be removed by coupling, then the 
Mathieu mode will grow very large as if an oscillatory 
load were being applied. 

Although the results are somewhat limited in number, 
they indicate that the response of a particular flexural 
mode has little direct effect on the response of another 
flexural mode; that is the intermodal coupling of the flex- 
ural modes is weak. However, all of the flexural modes 
have a very significant indirect influence on the response 
of the Mathieu modes, since they extract energy from the 
axisymmetric mode. Consequently the response of all the 
modes must be computed simultaneously in order to obtain 


the total response of the cylindrical shell. 


CHAPTER, Vi 
SUMMARY, FINDINGS, AND RECOMMENDATIONS FOR FUTURE WORK 
I. SUMMARY 


A numerical algorithm was developed for computing 
the dynamic response of a ring-stiffened, nearly circular 
cylindrical shell under transient axisymmetric loads. 

The governing equations, which are based on the nonlinear 
Donnell theory, account for the class of radial imperfec- 
tions whose circumferential variation can be expressed 

in a Fourier cosine series. The equations for the rings 
and shell structure were derived by using Hamilton's 
principle and include the radial inertia and a viscous 
damping term. The analysis accounts for the discrete 
location and eccentricity of a set of identical rings. 

Elimination of the circumferential variable was 
achieved by assuming a Fourier sine or cosine series for 
the displacements. The nonlinear terms generated products 
of Fourier series which were expanded to single series, 
yielding modal coupling terms. The series were truncated 
for the actual computations, with a maximum of fifteen 
participating modes. The retained modes were not limited 
to the first fifteen, but were selected to include those 
modes most responsive to the particular load. The partial 
differential equations were replaced with difference- 


differential equations by approximating the axial 
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derivatives with finite differences. A nonconventional 
differencing procedure, commonly called "modified" or 
"half station" differencing, was used. Provisions were 
made for both the clamped and simply supported end 
conditions. The iterative Newmark beta-method was employed 
for the time integration of the radial difference- 
differential equations. The axial and circumferential 
equilibrium equations were solved by Gauss elimination 
during each iteration. 

The algorithm was used to compute the response of 
a Simply supported shell under a uniform, exponentially 
decaying, radial pressure. This particular shell and simi- 
lar loading conditions have been considered in the earlier 
investigations of Anderson and Lindberg [19] and McIvor 
and Lovell [20]. The load parameters were selected so as 
to fall in the regime where dynamic buckling is governed 
by an elastic model, and the peak pressures and total 
impulses considered were between the static buckling limit 
and the dynamic buckling limit. Several response histories 
were computed, both with and without viscous damping; and 
comparisons were made with the results of References 
(19, 20]. The significance of the intermodal coupling of 
the flexural modes and the feedback coupling co the axi- 
symmetric mode was also investigated. Finally, the role 


played by the imperfection magnitudes was briefly considered. 
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A potential user of the computer program, or one 
interested in pursuing this work further, should be aware 
of the following limitations of the program: 

1. The analysis employs the Donnell shell theory 

which is not valid for many classes of shells 
[48]. 

2. The circumferential variation of the initial 
imperfections is expressed in a Fourier cosine 
series, instead of a complete series containing 
both sines and cosines. 

3. The axial, circumferential, and rotatory inertia 
of the rings and the shell are neglected. 

4. The portions of the computer program that cal- 
culate ring effects have been verified in only 
a cursory manner. 

5. A maximum of three rings can be accomodated by 
the computer program. The rings must be iden- 
tical in geometry and material properties and 
Spaced so as to divide the shell into segments 
of equal length. A principal axis of the ring 
cross sections must be parallel to the axis of 
the cylinder. 

6. The imperfection of the shell must not distort 
the rings. | 

7. The rotations of the rings must be small, and 
restrictive warping is not accounted for in the 


expression for the strain energy of the rings. 


ey 


8. The distance between the midsurface of the shell 
and the centroids of the ring cross sections 
must be small compared to the radius of the 
cylinder. 

9. The circumferential stress is computed only at 

O 


8 = 0 The meridional curvature change is 


not included in the calculation of the cir- 
cumferential stress; in addition, 1] -Vv : is 
taken equal to unity. 

10. In the axial equilibrium equation, the effect 
of a ring is felt at a station which is 


located one-half of the mesh spacing from the 


Gale. 
Lil —ebN DINGS 


The exponentially decaying pressure pulses were found 
to excite flexural modes belonging to two families, herein 
referred to as hyperbolic and Mathieu modes. The hyperbolic 
modes are the buckling modes, and they grow to their maximum 
amplitude in an exponential manner with little accompanying 
Oscillation. On the other hand, the Mathieu modes grow in 
an oscillatory manner requiring several cycles to attain 
their maximum amplitude. 

The hyperbolic modes grow primarily as a result of 
the quasi-static displacement of the axisymmetric mode. 

The loads considered in the present study were on the shell 


for a period of time much longer than the natural period 
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of the axisymmetric mode. Hence, the quasi-static response 
of the axisymmetric mode was nearly proportional to the 
instantaneous magnitude of the pressure. A given hyper- 
bolic mode was found to respond as long as the instan- 
taneous magnitude of the pressure exceeded the static 
buckling pressure of the mode. 

The Mathieu modes were found to be those modes whose 
natural frequencies were approximately one-half of the 
axisymmetric natural frequency and usually much higher 
than the frequencies of the hyperbolic modes. The Mathieu 
modes were identified through the use of a modified 
Mathieu stability chart. The modification replaced the 
natural frequency of the flexural mode by a pseudo fre- 
quency which accounted for the instantaneous value of 
the pressure. Hence the susceptibility of a given flex- 
ural mode to parametric excitation depended not only on 
the peak pressure and frequencies of the flexural and 
axisymmetric modes, but also on the instantaneous value 
of the pressure. Thus, as the pressure pulse decayed, the 
points corresponding to the flexural modes moved across 
the Mathieu diagram. Some modes moved from an unstable 
region into a stable region, while other modes moved from 
a stable region into an unstable region. The response of 
the Mathieu modes could be characterized as beat-type 
oscillations during which energy was cyclically exchanged 


between modes. 
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Anderson and Lindberg [19] investigated the dynamic 
pulse buckling of cylindrical shells subjected to axi- 
symmetric loadings. They showed that dynamic buckling 
occurs aS a result of the rapid growth of the hyperbolic 
modes, and that the Mathieu modes are of little importance. 

McIvor and Lovell [20] studied the nonlinear dynamic 
response of shells under purely impulSive loads. Their 
results show that for such loads only the Mathieu modes 
respond Since the quasi-static response of the axisymmetric 
mode is zero. However, the stresses developed during the 
response were found to be several times larger than 
stresses of an axisymmetric response to the same impuleuge 
load. 

The results of the present work consolidate the find- 
ings of these two previous investigations. At the lower 
peak pressures and total impulses, for which the quasi- 
static response of the axisymmetric mode was small, the 
Mathieu modes dominated; and the computed motion was quali- 
tatively similar to the motion computed by McIvor and 
Lovell. As the peak pressure and total impulse were 
increased, a gradual transition occurred; and the hyper- 
bolic modes played an increasingly significant role. Near 
the dynamic buckling threshold, the hyperbolic modes com- 
pletely dominated the Mathieu modes. The dynamic buckling 
loads predicted by the fully coupled nonlinear analysis 


agreed very closely with Anderson and Lindberg's results. 
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However, at peak pressures and total impulses well below 
the dynamic buckling threshold, stresses in excess of the 
(enc idee t eS ireSsmoccurredwduceeouthe icombined response 
of the Mathieu and hyperbolic modes. 

For those solutions in which a small amount of damping 
was assumed, the maximum amplitudes of both the dominant 
hyperbolic and Mathieu modes were reduced by 10 to 30 per 
cent in comparison to the zero damping case. Generally, 
the Mathieu modes were attenuated more than the hyperbolic 
modes. A critical pressure-impulse curve based on the 
yield point stress was determined for the damped response. 
This curve shows that stresses in excess of the yield point 
stress can occur at total impulses and peak pressures less 
than those required for dynamic buckling even when damping 
1s accounted for. 

The intermodal coupling of the flexural modes and the 
feedback coupling to the axisymmetric mode did not signif- 
icantly affect the response of the hyperbolic modes. 
Furthermore, the maximum amplitude attained by a given 
hyperbolic mode was found to be linearly proportional to 
Shape and magnitude of the imperfection in that mode. These’ 
findings are in complete agreement with the assumptions made 
by Anderson and Lindberg eon. 

However, the feedback coupling to the axisymmetric 
mode very Significantly affected the response of the Mathieu 
modes. This coupling is the mechanism by which energy is 


exchanged between the axisymmetric mode and the flexural 
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modes. On the other hand, the maximum amplitudes attained 
by the Mathieu modes were unaffected by both the inter- 
modal coupling of the flexural modes and the magnitude of 
the imperfections. 

The findings presented here apply to one specific 
Simply supported shell and loading condition. The buckling 


mode shape of this shell is sin—cos66, and the dominant 


L 
Mathieu mode is sinscos146. For other shells with diftipren— 


ent loading conditons, hyperbolic modes, and Mathieu modes, 
the coupling may be more significant; and these conclusions 
would not apply. This remains to be seen. 

The algorithm developed in this dissertation is 
capable of computing the fully coupled, nonlinear dynamic 
response of cylindrical shells under transient pressures. 
However, the practicality of the method is somewhat limited 
by the execution time of the computer program. For example, 
approximately forty-five minutes are required on the 
IBM 360 model 67 to compute a solution through thirty 
cycles of the axisymmetric mode when fifteen Fourier modes 
are retained. When the pressure pulse decays rather slowly, 
as it did in cases B and D of Chapter V, the solution must 
be computed through as many as fifty cycles of the axi- 
symmetric mode in order to encompass the response of all 


of the Mathieu modes. 
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III. RECOMMENDATIONS FOR FUTURE WORK 


The computer program has yet to be fully exploited; 
and, aS previously pointed out, the correctness of those 
portions of the program which compute the ring contributions 
have not been fully confirmed. Consequently, the following 
studies should be conducted using the program: 

1. Additional solutions carried out to longer times 
are needed to establish the significance of 
the intermodal coupling to the response of the 
Mathieu modes. If part of the coupling terms 
could be eliminated, the computing time could 
be reduced appreciably, thus, extending the 
usefulness of the program. 

2. Solutions should be generated for shorter shells 
with L/a = 1. For shells of this geometry, 
the static buckling mode and the Mathieu mode 
coincide. The dynamic buckling loads might 
turn out to be significantly different from 
those predicted by Reference [19]. 

3. Solutions should be computed for shells with 
clamped ends and compared with the results 
for the simply supported shell. 

4. The stabilizing or destabilizing effects of ring 
stiffeners should be studied. For example, if 
the stiffeners are too close together, it is 


possible that coincidence of the Mathieu modes 


di? 


and the hyperbolic modes may render the struc- 
ture unstable. 

5. The program should be used to study the dynamic 
Stability of shells under moving loads. 

6. An effort should be made to improve the Fortran 
coding of the program. It is believed that a 
professional programmer could improve the 


coding and, thus, reduce the running time. 
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APPENDIX A 


INCONSISTENCY OF THE CONVENTIONAL DIFFERENCES 


We will first show that the difference equations 
obtained by discretization of the governing differential 
equations do not agree with the difference equations ob- 
tained by minimizing the finite difference form of the 
potential energy when conventional central differences are 
used for approximating the derivatives. We will then show 
that the two sets of difference equations are identical 
when modified differences are used. Chuang and Veletsos 
[31] presented a very similar proof of this inconsistency 
of the conventional differences, but they did not expand 
the displacements in Fourier series as has been done here. 

Consider the static equilibrium of a cylindrical 
shell under a uniform radial pressure. For the purpose 
of the present discussion, the linear shell equations will 
be used. After deleting the acceleration and velocity 


terms, equations (3.15) become 


ee i ay nV, ; = = nu" - VW ¢ = 0 
ae - = nu, + =v" + nw = 0 
(A.1) 
ee 2 aaa 2 an, _ 
Ty © AW enece en a eee 


n n n 
Ve Weel 
S onP ae 
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The conventional central differences for the first, second, 


and ™rourth derivatives of a functions £(&) are 


- 


Te ame ore 
yee _ 
Wipe 3 (cece He ) Cee) 
ee a a" 
‘Preece x ~ 5F (xta 7 “qtr * OfK 7 AE * FA 2) 


where d is spacing of the node points. 
Using the conventional central differences given by 


(A.2) to approximate the derivatives in equations (A.1) 


leads to 
3 pe ee eg MeL) - apr 8 Us 
oa (Weed Weep) ie 
-n°vp = = n : (ie, - Ue = = We 20 KV) nw, = 0 
= a" I=y (we -4we,, + Owe - 4we_, + wes) (A.3) 
on 3 (wes - 2w.. + Wey) + now? = 


A set of difference equations can also be derived by 
minimizing the potential energy. The potential energy II, 


for the linear shell under consideration here, is given by 


iS 


L 2%a 


a 
IT == N oy + N 
a Ey 
(A.4) 
= MMe + OM kay = sass =) 2ZPWidyax 


The strain-displacement relations of the linear theory are 


2 = ae ete ee W, 
xX xx 
ge =V,_- uM x = W, (Ans) 
MY ry! a ye ey 
Y xy = Ury + ae Xsey = Wry 


The in-plane forces and moments are related to the 


displacements by 


_wW 

Ny = B[U,, - ie =) 
N = B(V,. - — + vU,_) 

Y Y x 

1l-v 
Go ae 
(A.6) 

M = -~D(W, vi 

xX xX yy 
M = -D(W, + av, ) 

Y YY xX 
M = (l-v) DW, 

XY xY 


As before, the displacements are expanded in Fourier series 
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Wisc) = -£ U(x) cos né@ (A.7a) 
n=0 

Wiser) = oD Vv (x) sin né (A. 7b) 
n=l 

W(x,8) = 2 We (x) cos nég (A.7c) 
n= 0 


The terms comprising the integrand of (A.4) can be 
expressed in terms of the displacements by means of (A.5) 


and (A.6), leading to 


N € = B[U,- taviiVv,.- wl io, | 
x x a x 
- W 
N € = BL(V - =- Te meals igs 
yey ME OU aa ale Ong 
l= 2 
xy ' xy zB (U, a V, ) 
(A.8) 

ioe = ae (Ww? + vW, ’ ) 

ox xX ce VV 

= 2 

M x a -D (W, vW, W, ) 

ney, yay, YY xXx 
oe PST oy Dw 

xy XY "xy 

Replacing the displacements in (A.8) by the series of 
(A.7) yields 

ieee = 2 Una es ny = wu, 

mx x Xx a x 

n=0 m=0 
cos n@ cos m6 (A.9a) 


137 


oO 


Ne, = BE © [+ (nv - w)(mv™ - wy + 
yY*  n=0 m=0 a 





< Ure (mv ~ Ww") ] cos n8 cos m@ (A.9b) 
law . a ny: Jee mi _m™ .m, , 
Nyy’ xy 2 ees ea ee a a ie a J f 
Sin n@ sin mé (A.9c) 
Myoe spr or wel? - See 
x xX s 7 xx Kx Z xX 
n=0 m=0 a 
cos n8 cos mé (A.9d) 
= = n“m@ wwe a nym 
MX, = “Dz 2 uns - v5 WW,..)° 
YY n=0 m=0 a as a 
cos n8@ cos mé@ (A.9e) 
Me = (10) ene 
xy “Xy Aa ay 42 tek 
Sin n8 sin mé@ (A.9f) 


The potential energy of each mode is obtained by substitut- 
ing (A.7c) and (A.9) into (A.4) and performing the circum= 
ferential jntegration. When m, n # 0, the potential energy 


tela 


of the n mode reduces to 
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2 x 
Vv n n n Vv 2 
SE Aca Vieeeea juliet DI (W,-) 2 a I 
n° - 
+ 2(l-v) D — (W,_) (A.10) 
a 
n n 2 na eee 
+ Dia (Wo) - v5 WW, ra dx 
a a 
When m = n = O, the potential energy is 
L 2 Zz 
Bes 0 eee CO i. 0 
Ml = ape olh {B[(U,,) = W oe a" > (Ww) J 
0 a 
0 2 0 
aD ir eee mea cg } ax (A 11) 


The total potential energy is given by 


The principle of Minimum Potential Energy applies to each 
mode individually as well as to the collective sum of the 
potential energies of all the modes. Furthermore, the 
potential energy of a mode may be multiplied by an arbitrary 


constant without changing the conditions necessary for its 
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minimization,, “Accordingly, n° is multiplied by one-half. 
This allows us to write a general expression for the 
potential energy which is valid for all values of n. For 


convenience the terms are grouped into two integrals, 


leading to 
I = may {B[ (Uy )" + 2 (oy wo 
2 0 ’x a ox 
ne wy 4 + D[(W,™ )* =o ee nee 
2 ’xXX Z SCX 
a a 
te wn 7 n 
* oa (Wo) J + 255 ,PW } ax (A.12) 
a 
L 2 
Ta 1l-v n ie a 
ee 
0 
2 2 
n n 
+ 2(l-v) D > (W,,) } ax 
a 
where Son = 1] when n = O and Son = 0 when n # O. 
Let I be defined by 
=n 2 n 
oe —<—s II (Ay bsp 


Expressing (A.12) in terms of the nondimensional 


variables of (2.21) and using (A.13) give 


a U 2 
T= fs { (ure) + 2vi(nv" = wit) eer (nv™ = wh)? 
0 S 
DT oe a eon ne ar ene os pw} dé 
12 BEE "ED On 
‘ ; (A. 14) 
+ ; =. [vee _ nu’) J + z (1-v) a?n® (we) \emene 
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Conventional differences will now be used in sengune: 
tion with equation (A.14) to derive a set of difference 
equations. In order to discretise the potential energy 
integral, let the € - coordinate of the a mode be divided 
into M elements each of length d. At the center of each 
element, locate a node point K. Then ie may be expressed 


by the following series. 


M 2 
=n n n n ie 
t= Dimes (Ups) + 20  — we) (a, ) 
Mo K=1 = K i 7 a 
2 2 
n n 1 2 n 2 an n 
+ (nv, - Wy) + I> IWreedy - 2vn We lWiee)y 
Aare ee n 1l-v n 
+n (w..) ) + 25 5 .PRWK + a U(vigdy - (A.15) 
ae tev oe ee} a 
K 6 "EK 


where d = 1/M. The potential energy I can be approximated 
by taking M to be finite and using central differences to 
express the derivatives. 


The principle of minimum potential energy demands 


that the displacements ee, Vie wae ue ey Ee eae ae Wa 
minimize the potential energy Il’, i.e. 
M =n on =) 
sr = on bux + on évy + on bw = 0 (A.16) 
K=l1 du dV dW 
K K K 
. ’ ; n n n 
Since the virtual displacements buy 6V7 1 OWy ree OU OViye OWay 


are arbitrary, it follows that 


on 


= 





Q 
a 
AS 


a 





— (A.17) 


@ 
< 
AYP 


- 





Q 
a 
AD 


When the derivatives in (A.15) are replaced by their con- 
ventional central difference approximations (A.2), the 


potential energy becomes 


n 


u - rs Z uo wu? 2 
—n K+2 K K K-2 
= ooo 4 a . 2 se 
I { ( 54 ) ( 4 ) 
UKE? ya 
n n n a 
t2v (nv, 4 - ee + 2v (nv, - Wy) 
n n Nun 
pane ERD) Py inv? et ( 5) 
2a K-1 K-1 2d 
n a n n 
a aoe 2 Kt2 viet ee ae 
K K LY. ge 
n _ n n Tae n n 
Wied Zw, + Wee} 2 Wy 2Wyiy + Wee 2 
+(e) + on 
a? 2 
n n n rw 5, n 
2 on Weta ~ *¥x4 1 + WE 2n “Rel °K? “eed 
-2vn~ w ( 2 eee ) 
K+1 2 K Z 
n n : n : 
WwW. - 2W + w 
2n K K-l K-2 eee Tietz 
-2vn Wei} ( BZ )teeetn (wi) )+-° +205 WPRWK 
no = yn n n- 
grey (Kt? K nut 5? 4 (Kt YK=1 ae 
2 2d K+1 2a K 
voy a we we D 
K-2 2 2 K+2 K 
a — +—(]- tooet 
( 4 nuy_,) q (1 ee reel 4 ) 
we we 
~ ee 
(A) “1 + } 4 (A.18) 


eZ 


Only the terms containing es Vi and We have been displayed 

Since the displacements at the other stations, Kt+tl, K-l, 

etc., do not enter into the variation of Ne, Te and Wie 

The axial, circumferential, and radial difference equations 

are obtained by equating to zero — a, an The result- 
by) wd hy 


ing equations are 





Axial: 
mn - 2u” + u? ve - v2 
vee K o +V + J 
ns sca * n( Ae) iouaioen 
4d IMG | 
we W 
=) a 
-+ n-yn a Rae K 1) <9 
Zea 
Circumferential: 
u® u? ve Dice 
2eneylAVeeKe DMeK=liyedlinwey) Ke2 Me Kee, 
s 7 2a 2 4q2 
2 aa 
+ nw = 0 (A.19b) 
Radial: 
n n n n n 
eer > a ew aw Pw 
TD a [{ 7 ) 
a 
n n n 
W - 2w.. + w 
on? | K+2 K = oe wh i 
(A.19c) 
n n n 
v2 ceo Sey te Ob ~ Wy gee > 
> nd ( A 0 
a 
un u” 
6 nN ae = K+1 K-1 
ee ee ae 2d 
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Let us compare the above difference equations with 
(A.3) which were derived directly from the differential 
equations. Note that the second derivatives in equations 


(A.19) are represented by expressions of the form 


1 


EoD Elm Ke 


2) eee 
4d “ 


Ra: 
aS opposed to the conventional central differences of the 
second derivatives in (A.3). 

In addition, (A.19c) contains a term which has no 


counterpart in (A.3c) namely 


n n nN A) 
vy (2.2.2 (“Keo “Wks * °WK ” ““xeamalell 


am ( 


Were it not for the pressence of ae this term might be 
discarded as a higher order term. However, the problem 
under investigation encompasses response modes having 
several circumferential waves; and the term may not be 
small compared to other quantities which appear in the 
radial equilibrium equation. 

Consider now the modified central differences. A 
set of difference equations will be derived by minimizing 
the potential energy in modified difference form. Conse- 
quently, a series representation of (A.14) must replace 
the integral formulation. The right side of (A.14) was 
written as the sum of two integrals for reasons which will 


now become apparent. In order to obtain a series which is 
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equivalent to the first integral, let the €-coordinate be 
subdivided such that a whole station node is located at 

the center of each segment. This procedure is schematic- 
ally illustrated in Figure A.la. In approximating the 
first integral, the half station at each end is neglected. 
Doing so does not affect the difference equations that must 
hold at the interior points. Hence, the first integral 


san be written in the form 


M-1 2 2 
. n n> <n n me 0 
lst Integral = 7 {(are) y + 2v (nv, Wy) Ure) y + (nv, Wee) 
Pes? hie ee Sune wl oa De) ee Deepen a 
12 SEE Ok K VEE ak K On? K K 


The second integral of (A.14) is approximated by using 
finite elements which are centered at half station nodes, 
as shown in Figure A.lb. The series approximation of this 
integral 1s given by 


M-3 


2nd Integral = (Sti, 2g - aul? 
K= 1/2537 2». 28 
ja Gls, see ee ye d 
6 EB) K 


Combining these two series leads to the equation for the 
potential energy. Replacing the derivatives by modified 
central differences of (3.16) and minimizing the potential 


energy lead to the governing difference equations 


os 


FIGURE A.la 
“WHOLE" STATION MESH-MODIFIED DIFFERENCES 


M-I/2 





ELEMENT +F 2 


FIGURE A.Ib 
"HALF" STATION MESH- MODIFIED DIFFERENCES 


PG 


n n l-v hk n n 
_ on = ou = 0 ) +——n = (v - V_,) 
—s + K 
42 (UK+3/2 K+1/2 K+1/2 2 d K+1 
a) a Peet) ie 
ome i emma Kate Ko 
eee. 1 1l-v An ere ea Pi, n n 
ei ee coum > | oe eo Kl OK ene 
+ nw! = 0 (A.20) 
K 
1 ee n nN n rn 
di ace eK i i) 
2 i n n n ag 
2n = ee a Yee on ae | 
S i. ee digger naa 
Cnr ee ee iy Ken 2) 


These equations are identical to the static version 
of (3.17), equations obtained by applying the modified 
differences to the governing differential equations (3.15). 

When the modified differences are used to approxi- 
mate the derivatives, the set of difference equations 
obtained by satisfying the differential equations at a 
finite number of points are the same as the difference 
equations obtained by satisfying the condition of minimum 
potential energy. On the other hand, when conventional 
central differences are used, the difference equations 
obtained from the differential equations do not agree with 


those obtained by the energy method. Consequently, the 


Heo? 


modified differences should yield more accurate solutions 


than conventional differences for the same mesh length. 
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APPENDIX B 


ADAPTATION OF MODIFIED DIFFERENCES 


TO THE NONLINEAR EQUATIONS 
I. EVALUATION OF NONLINEAR TERMS 
In-plane Forces, )\1, and A2 


The in-plane forces were expressed in Chapter III 


bye(s.5). Let ne and ng be defined as whole station 


quantities, and let n? 


Ee 


Replacing the derivatives in (3.5) with the appropriate 


be a half statwven quantity. 


modified difference expressions from (3.18) leads to 


a 
(np) x = (uy, - ae) + vinvy - Wie) 
(Bea 
n n n n 
ae ou t+ MC ae ) 
XX), itl. TT. Itt. 
7 ees v eee art 
(Ne), = nv We + 5 (uy U,_4) 
(B. 1b) 
n n n n 
+ Cc + C + v(C + C ) 
Try ITT. AX, IXX,, 
n lve. F n Paes n 
leak oo hae Ke et om 
- ee (B.1lc) 
+ C + C ] 
TXT) ITX 


The nonlinear terms in the axial and circumferential 
equations are contained in hl” and ae defined by (3.8) 


adiGwelo.oe  Hewatton) (3.20) requires that 11” be defined 


Eo: 


at half stations, and equation (3.21) requires that 
h2" be defined at whole stations. ADProOx<tMaet ing wears) 
at the half station k and (3-9) wat the whole Stavpiene 


TesWlecunm 


ele n lap n 
hl, = =<e = Cl) + =e C eee 
Ie d XX ay XX, d IXX, 41 IXX, 
1 sin n i el n 
v{=(C —se= je =e =-C 7] (B.2a) 
G TT. 1 Ty. d ITT 4d Ty. 
1-v n n n are Tl n 
fee Ce | a ee 
k ire k 
ogc n n n n 
h2,, n [Crm + Crop tv (Cy +Cryy )j+ 
K K K 
l-v ,l n 1 n hn 
—— [= (C -C )+ ce =-C ) (Baza 
2 ad XxX, "do ee oe 
5 = (Ce 7 ae ) nWy 
k k-1 
; : n n 
Examination of “%Bll) and (68> Z2jeshevs thae Cras Crimp 
n n ; ; n 
Cyys and Cryy must be defined at whole stations while Cryms 
nN nN « 6 
Cyms and Cimy must be defined at half stations. As seen 
n n ; 
from (3.4a) and (3.4d), Coy and Crux are obtained by 


summing appropriate products of the form 


and 
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n 
IXX 


and we must be calculated at whole stations. However, w 


E 


is defined at whole stations only; and the modified 


Since eee and C are defined at a whole station, Wie 


difference form defines Wire at half stations, not whole 
stations. A similar situation exists for the half station 
nonlinear terms am and oo terms which can be evaluated 
only if w is known at half stations. This creates some 
difficulty when utilizing the modified finite differences 
to solve the nonlinear equations numerically. The diffi- 
culty 1S overcome through the use of interpolation formulas 
or simple averaging. 

Consider the problem of determining Wie at a whole 
station. With reference to Figure B.1, let Weed and we 
be the radial displacements at half stations k-1l and k, 
respectively. The required derivative of as which we 


shall denote (w', 1s given by 


ghei 


(Ee) 


By expanding w in Taylor series, the quantity (wy = we a) 


can be expressed as a linear combination of Wee We yr 
n n . : n 
Wea! and Wea: Using the Taylor series to express Wir we 
have 
ww Siwy 4 Bye) + 0 (83) (B.4) 
ie ie 2 ame k 22 3 


Zee 


a W + 
k-1 k 
k-1 k 
K-2 K-1 K K+1 K+2 
w rr n aot Ee 
m7 K-1 ; K+1 K+2 
FIGURE B.1 


FINITE DIFFERENCE STATIONS 


n 
where (w, 


and (w,,-), are the first and second deriva- 


n 
bade Eb )e 
tives of w’ at half station k. The notation 0 (a>) means 


that the truncated terms are of order a>. The second 


derivative (Wire), can be expressed in terms of (Wrep)y OF 
ree) eed a 
n 2 n den 2 
(Wreedy = (Wreedy it a Wreeedy + Ona) (BStay 
(re) pa Wie egy > Siw eee (B. 5b) 


Adding the above two equations yields 


wii = sliie * Wit er wae 
5 5 Kee - Wear 7 WK + MK) + 0(a*) (B.6) 
Also, we have that 
wre) = Bey 7) + 08) (B.7) 
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Substituting (B.6) and (B.7) into (B.4) leads to 


Nees 2 
W SS = 


n n ah n 3 
K = eK + Mea) 7 TER + Meg) + Old) (B.8) 


It follows immediately from (B.8) that 


n eet n 3 
+ wy) - 7¢elWeiy + Wy_5) + O(d eGo) 
Finally, the interpolation formula for the derivative is 
obtained by substituting (B.8) and (B.9) into (B.3), 


resulting in 


n el Dee ee 1, - eek 2 
(Wi pd = Teq (1° (Wray Wy) (Weis Wy_>)] + O(d) (B.10a) 


nat a half station, 


When it 1S necessary to calculate w 
either (B.8) or (B.9) may be used. 
At the boundaries the first derivative of w” may be 


approximated by using a forward or backward difference with 


(2.8) omm(B.9)@e Since wy = Wa = 0, the result is 
n eee ee n 
(Wredyoy = (a) lye Wo iG + wa] + 0(d) (B.10b) 
and 
te . . Acari lee n 
(Woe) von = (3) (Fe ¥u-1 eee wy) ena) (Bade) 
The procedures for evaluating the Fourier coefficients 
of the nonlinear terms defined by equations (3.4a) - (3.49) 
are described in Appendix B of Reference [30]. For example, 


Cor at whole station K is given by 
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n Jieeemrediy c, itn C i-n| 
= = + 
(Cyd = | eer rom Wien, MmCe Rite Ke] (Be11) 
where 
Pa ete 0 
é 
NY — 
0 if n= 0 
tit 1 Aen 
Cc 
U = 


23ifie=zn 


The terms of the series can be evaluated by utilizing 


(3.18), (B.8), (B.9), or (B.10). The other nonlinear terms 
n n n n n n 
Cry ; Co . Crop ; Cyr ; Cryp , and Crrx can be evaluated 


K K K K K 1 K 


from equations similar to (B.1l1). 


Nn Ti 


n 
Nonlinear Terms "é Neg , and Ng 


The nonlinear terms in the radial equation, (3.22), 


are embodied in the terms n” . no ; ne ,and nn - Equation 
OK SK 58% OK 
(B.lb) 1S used to compute the first of these terms. The non- 


linear terms ne and ne are defined by (3.4h) and (3.4j). 
The calculation of (np) involves the summation of products 
K 


of the form 


i j = 3 
(ne) ,mliw? op em (w? 


1 n n n n 
XX," =i Toe Crypr and Crny are 


not truly nonlinear since they are Sums of products that 
contain derivatives of the initial imperfection. However 
they are computed in the same manner as the nonlinear 
Geminis 


Themcock te. crent smc 
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for appropriate values of i and j, and the calculation of 
(ng) requires summing products of the form 
K 


(nb) t-37w) + W5) 

K K 

The summations are defined by expressions similar to 
(B.11). The terms in the products can be computed from 
(B.la), (B.lc), and the difference approximations of 


(on Loe. 


The remaining nonlinear term (nie g) requires the use 
K 


of a simple average. Equation (3.41) expresses (neg) as 
K 
oO I ; OO 7 = ©O n 
(22 Sin 16)[ © -j (w;,-+w;7,)sinjOJ= =n cos né 
i=l 5x j=1 Sk n=0 5x 
(3.41) 


However, (B.lc) defines Neg as a half station quantity. We 
n 

Eo, 
can be written as 


ry 
a 


nN 


take n as the average of n and n * thus, (3.41) 


Seed 


cos n6é = [ £ (n + nt \Sinmeo} |. > =e ap ) 
“=1. 68, 8 Oy j=1 Ex §&x 


Sin 36] (B. 12) 


n s e e e Nn 
where, as before, w, 1s the first derivative of wat 


eK 


whole station K and can be evaluated by using equation 


(5. LO) se Eeom (Bale) itetol lows that 


n n =. 
n ton = i=v [ = 


Nn 
= ais 2 K+ K-1 eee! 


n n n n n 
+ C + Cee tC + C +C } 
IXT, _, “ITX, “ITX,_y 
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Finally (neg) can be evaluated by using a formula from 
K 


Appendix B of Reference [30], similar to (B.11l). 
II. Boundary Conditions on Displacements 


Consider the case in which the axial displacements 
are zero at the ends. Since u’ is not defined at the whole 
Stations K=l and K=M at the ends of the shell, interpolation 
formulae are required. The finite difference stations near 


the end € = 0 are shown in Figure B.2. The displacement 


Ue can be expressed by the Taylor series 

n n ayaa 1,d,2,0.n 2 

= -_ — —— ie of 2 

aa aaa 2, 5 ‘5? oP ae 0 (d ) (B 14) 

Now 
aC - ur 
uy = el gpm 2 + 0 (a?) (Beel!5 ) 
ae d 


and using forward differences gives 


n (Ure) yee = (Ure) yay 
(Ure) = ae ay ei + 0(d) (B.16) 
K=] 
We also have 
u” = ce 
(arp) = a Pod 4 | (B.17) 
K=2 


To obtain the desired boundary condition, (B.15) and (B.17) 
are substituted into (B.16), (B.15) and (B.16) are sub- 


stituted into (B.14), and (u”) is equated to zero. The 


K=] 
result is 
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"0 ma oa 
a 4 —_o—_+—_] 
K=] K=2 
FIGURE B.2 
FINITE DIFFERENCE STATIONS 
AT THE END € = 0 
n n 1 on Se. 
une 2u, 5 uy, + O(d~) = 0 (B.18) 
similarly, the condieron tier vu’ = 0 at —€ = 1 can be ex- 
pressed by 
n n ilieeaaela Bo e-= 
ut 2uny 7 z Uy» + O(d-) = 0 (B.19) 


When the shell is clamped at K=l and K=M, the slopes, 
(Wie) pay and (Whe) yom must be zero. Equating (B.10b) and 


(B.10c) to zero yields 


n n Hite. 

Wo = Iw. +. W = (0, (B.20) 
n nN n _ : 

Ses eee = Oo (B.21) 
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In equation 


APPENDIX C 


DEFINITIONS OF CONSTANTS AND MATRICES 


Core.) 


In equation (3.23) 


= nc 


n 2 
Ae [C,e (2-en ) = Ce] 


een’ - c, (n*-1)? _ C, (en*- 1)* - 2d 


4 6 


2 
n 


ae [C, (l-en*) + Cel 


d 


3 (l-en*) 


In equation (4.2) 
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ii dn (—— 
oe am 
A” = 
l-v 
s 2 
Pee Nr pep 2 
-V 
= 2 da“n om) Z Say. dén'(C,n°+C,) ~ (2%) an 
j=1 3 2 
re 
By = 
2 Ne 2 
- aay dn -l(l-v) - a n* + Zz 6 a°n Cc 
2 eas 3 
J=1 
1 6) 
Ge = 
dn (==>) = 
J 


In equation (4.8), a is given by the following expressions. 


Simply supported ends: 


l+v 
1, Fug eee), ine aoe 
2 2 
n — 
Py = 


Clamped ends: 


L+v 
-4/3 id 
cary te ant 
n —_— 
Py = 
0 0 
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n 


In equation (4.11) matrices H and 0” are given by 


the following expressions. 


Simply supported ends: 


1g: d*n* (SY) dn (=5) 
H = - 
an (=) 1 sae 
1 0 
Qo" = 
eg = 
d(Cyy + Chyy) 
ee = Sha | + 
0 
M 
Clamped ends: 
ha nee any) 
2 2 
Ho = - 
an (=5%) han aden 
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Oe = 
1+v 1-v 
jan a, 
nen 
a eed 
In eguation (4.12) 
0 0 
R= 
1l+v l-v 
dn (=) ae 
1l-v i gen? (=a Ht gen? 
2 2 
0 
n 
S ans 
me 
m 
a¥ in azn? 


Ze 


APPENDIX D 


COMPUTER PROGRAM 


I. FORTRAN VARIABLES 


Table XI lists the major Fortran variables and their 
algebraic counterparts from the main body of the disser- 
tation, if such a counterpart exists. The equation in 
which the algebraic counterpart was first defined is indi- 
cated in the remarks column. The more important program 


control variables are also briefly described. 


Il. ‘DESCRIPTION OP SHE PROGRAM 


The computer program has been written in FORTRAN IV 
language for operation with the IBM 360 MODEL67 Computer. 
The program was compiled by the FORTRAN IV G Level l, 

Mod 3 compiler. The program is dimensioned to accomodate 
twenty equally spaced axial stations and fifteen arbitrary 
Fourier modes. The program will handle up to three 
indentical ring stiffeners, located so as to divide the 
shell into segments of equal length. The memory require- 
ments and the execution time are contained in part IV of 
this appendix. 

The program solves a set of three nondimensional, 
nonlinear equations by an iterative procedure. The depéens 
dent variables are expanded in Fourier sine or cosine 


series in the circumferential coordinate. The nonlinear 
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terms are thus converted to products of Fourier series. 
The products are expanded to obtain a single series for 
each such product. Derivatives with respect to the axial 
coordinate are evaluated numerically by uSing a set of 
modified staggered node, central differences. 

The radial equation of motion for each mode is inte- 
grated in the time domain by using Newmark's beta-method. 
After determining the radial displacements, the nonlinear 
terms in the axial and circumferential equations are 
evaluated. The axial and circumferential displacements 
are then computed by Gauss elimination. Finally, the 
radial accelerations are computed from the radial equilib- 
rium equations; and a convergence test is made. The 
process iS repeated until convergence is achieved or until 
a specified maximum number of iterations have been com- 
pleted. If the solution fails to converge in the maximum 
allowable number of iterations, the time increment is 
multiplied by two-thirds and the procedure is repeated. 
Whenever convergence occurs in less than some specified 
minimum number of iterations, the time increment is in- 
creased by 15 per cent. 

The output consists of the maximum circumferential 
stress at 0 = 0° attained during the solution together 
with the time and station at which the maximum stress 
occurred. In addition, the displacements and in-plane 


Stress resultants are printed at specified intervals. 


Ze 


Main Program 
The main program performs the following functions: 
1. Controls the logical connections between the 
subroutines. 
2. Adjusts the length of the time increment. 
3. Sums the Fourier coefficients in the series for 


the in-plane stress resultants. 


4. Computes the circumferential stress due to 
bending. 
5. Computes the maximum circumferential stress 


together with the location and time of its occurrence. 
6. Prints out the solution at specified time 


steps and after completion of the solution. 


Subroutine START 

This subroutine performs the following functions: 

l. Reads in the data cards. 

2. Prints the problem description and program con- 
trol parameters. 

3. Sets the initial values of all array elements 
eo Zero. 

4. Calls subroutines WIDER 1 and WIDER 2, or 
WI1SIN and WI2SIN, which compute and store derivatives of 
the initial imperfection. 

5. Sets certain program control parameters to proper 


ihe | ese 


we 3 


Subroutine WDEF 
This subroutine performs the time integration by 


Newmark's beta-method. 


Subroutine MODES 

Subroutine MODES was taken directly from Reference 
[30] with the permission of its author. 

In MODES, the arrays that define those sets of indices 
that combine to equal each value of n in the problem are 
determined. They are IS, JS, ID, JD, IJS, MAXS, MAXD, and 
MAXSY. MODES is called during the first iteration of the 
first time step and during each subsequent iteration until 
the number of Fourier terms reaches a specified number. 

The arrays computed here are used as control parameters 
by subroutines LAMIN2 and TNETA in evaluating the nonlinear 


terms. 


Subroutine LAMIN2 

The quantities LAM1 and LAM2, the terms which com- 
prise the load vector in the Gauss elimination, are computed 
in Subroutine LAMIN2. The derivatives of the radial 
displacement that appear in LAM1 and LAM2 are first 
evaluated. A multiplying and summing procedure is then 
carried out to evaluate GxXx, CiIxXxX;, CLT, CilT, Cxt er 
and CITX, the Fourier coefficients of the expanded products 
of series. Arrays IS, JS, ID, JD, IJS, MAXS, MAXD, and 
MAXSY, which were generated in MODES, are used as control 


parameters in the multiplications and summations. 
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Subroutine PMATRX 

PMATRX serves as a control program for the computation 
and storage of the matrices P, DEE, and DST, matrices used 
in the Gauss elimination solution for the displacements 


u and v2. The subroutine is called by the main program 


k K 

during the first iteration and whenever the expansion of 
the products of series introduces new Fourier modes. The 
boundary conditions at the first axial station are in- 
corporated into the submatrices of P associated with the 
first station. Subroutines ABC and PANDD are then called 


to compute the rest of the P matrix and the matrices DEE 


and DST. 


Subroutine ABC 

This subroutine computes the elements of the A, B, 
and C matrices for a given Fourier number. These matrices 
are used in subroutine PANDD to compute the elements of 
P, DEE, and DST. At stations having rings, the matrix BR 


is used in place of B. 


Subroutine PANDD 
The matrices P, DEE, and DST are computed by PANDD. 
Subroutine MATMAT and INVMAT are called to perform ele- 


mentary matrix operations. 


Subroutine UVDEF 
Subroutine UVDEF computes the displacements uy and 


1) 
V 


KS which correspond to the new radial displacement, wl. 
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UVDEF first computes the new vector GEE, then performs 
row operations on the vector GEE to generate the X vectors, 


n n 
ke and Ve 


routines MATMAT, INVMAT, and MATVEC are called for ele- 


and finally computes u by back substitution. Sub- 


mentary matrix operations. 


Subroutine TNETA 

This subroutine performs the following functions. 

l. Computes TX, TXY, and TY, the Fourier cCoeffictenc. 
of the in-plane stress resultants. 

2. Computes WWX, WWXTH, and WWTH, the Fourier co- 
efficients of the &-curvature, the twist, and the circum- 
ferential curvature, respectively. 

3. Carries out a multiplication and summation pro- 
cedure to compute ETAX, ETATH, and ETAXTH, the Fourier 
coefficients of the nonlinear terms in the radial equilib- 
rium equations. The arrays IS, JS, ID, JD, IJS, MAXS, 
MAXD, and MAXSY, which were prepared in Subroutine MODES, 
are used as control parameters in the multiplication and 


Summation process. 


Subroutine ACELR8 

This subroutine solves the radial equilibrium equation 
of each mode for the radial acceleration WDDOT1. The con- 
vergence test is performed as the new accelerations are 
computed. The parameter NCONV is set to zero when the solu- 


ELOn converges, or to one, if the convergence test fails. 
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The program stops checking for convergence as soon as a 


tes@ fails. 


Subroutines INVMAT, MATMAT, and MATVEC 

These subroutines perform elementary matrix opera- 
tions on 2 x 2 matrices and 2 x 1 vectors. INVMAT inverts 
a matrix, MATMAT forms the product of two matrices, and 


MATVEC multiplies a vector by a matrix. 


Subroutines WIDER 1 and WIDER 2 

These subroutines compute and store the derivatives 
of the initial imperfection required in subroutines 
LAMIN2 and TNETA. The subroutines are based on the assump- 


tions that the imperfection in mode n is equal to 
wi {l- cos [2m (N + 1)&/i]} cos né 


and that the rings are located such that the shell is 
divided into ee 1 segments of equal length. This 

means that the imperfection and its first derivatives 

are zero at the ring stations and do not produce any 
strains in the rings, as assumed in the analysis. Accord- 
ingly, there is a finite number of compatible ring sta- 
tions. They are tabulated in the program listing of 
subroutine WIDER1. The program listed calls these 
routines. If another form of the imperfection is desired, 


a separate routine must be written. 


Subroutines WILSIN and WI2SIN 

These subroutines are substitutes for WIDERI1 and 
WIDER2 and were used in computing the solutions pre- 
sented in Chapter V. These subroutines are based on the 
assumption that the initial imperfection in mode n is 


given by 
Ties Sin(ne/i)=ecos. ng 


Subroutine PRLOAD 
This subroutine computes PLOAD at each station for 


a uniform time decaying external pressure, i.e. 
PLOAD (K) = PO exp (-RT) 


The subroutine may be replaced by the user for other 
pressure time histories. The variables PO, R, S, and 
OMEGA, in labeled common block Al, are read from data 
card number ten. These, variables are not used anywhere 
else in the program and are available to the user in 


developing an alternate subroutine for the loading. 


Subroutine AMPMAX 

Subroutine AMPMAX is not an integral part of the 
program. It was, however, employed in computing the 
solutions presented in Chapter V. It should be used only 
in conjunction with WI1SIN and WI2SIN. AMPMAX computes 
the amplification of those Fourier modes for which an 


initial imperfection was specified. If the magnitude of 
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the initial imperfection for some mode is set identically 
to zero by the input data, the subroutine will compute 

an amplification of zero even though the response of that 
mode may not be zero. The amplification is computed at 
the midspan of the shell or at the whole station nearest 
to the modspan. In computing the amplification, the 


initial imperfection is assumed to be equal to 


nh 
wi Sin ( cos nQ 


TE 
=e 


Iie. INPUT REQUIREMENTS 


The layout of each data card is shown in Figure 
D.1l. The first three cards will be referred to as 
Tite, 2, 3; and the ramaining caBuds,,as DAPA™1] - 913. 
The program uses nondimensional variables throughout 
and the user must supply the data in nondimensional form. 


Themcontents of each Icard are Iisted below. 


Data ards 
i Tee = oe The description of the problem is placed 
on these three cards and will be printed 
at Gop o& the Mirst page of output. 
DATA 1: 
MM - The number of axial stations, including 
the stations at the ends of the shell; 


MM @ 20 
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NR 


INMAX 


MAXM 


DATA 2: 


NITMAX 


NITMIN 


MAXSTP 


RoR INT 


TFINAL 


BETA 


Itesmunber Of Stiffening rings; NR< 3 
The number of modes with an initial im- 
perfection, including the axisymmetric 
mode (n = 0); INMAX < 15 

The maximum number of Fourier modes to 


be considered; INMAX < MAXM < 15 


The maximum number of iterations per time 
step (recommended value = 4) 

The minimum number of iterations per time 
step (Recommended value = 3) 

The maximum number of time steps to be 
computed 

Time step intervals at which the dis- 
placements, stress resultants, and maximum 
stress are to be printed. If KPRINT = 50, 
print out will occur at the end of time 
Stepus0; LOO, Ps077etc. 

The nondimensional time at which the 
calculation is to be terminated. (The 
nondimensional period of the axisymmetric 
mode is approximately 27) 

Integration parameter of Newmark's beta- 


method (Recommended value - .1666667) 
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EES 


DATA 3 


NBC 


NSYM 


DATA 4: 


- Convergence téesStppaarameter., § Convemaence 


1s assumed when the magnitude of the 
change in the acceleration for two 
successive iterations is less than or 
equal to EPS times the magnitude of the 
larger acceleration. (Recommended value - 


- 0100 ) 


The boundary condition parameter. When the 


Shell: 
1. Has simple boundary conditions, 
NBC = l. 
2. Has clamped boundary conditions, 
NBC = 2. 

Symmetry parameter. When the motion: 
1. Is assumed to be symmetric with 


respect to the midspan, NSYM = l. 


2. Otherwise, NSYM = 2. 


This data card contains nondimensional 
parameters representing the stiffness 
and geometric properties of the stiffen- 
ing rings. This card must be omitted if 


the shell does not have stiffening rings. 
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el 


CZ 


es 


C4 


G5 


DATA 5: 


DATA 6, 
DATA 7, 
and 


PATA Ss: 


I (1-v*) 


ach AX 


tly 


ae 
hAxX 
Z 
A (1-v ) 
hAxX 
L (1-v*) 


a“h AX 


2 
E Cp (i-v ) 


Qa°h AX 








where 
x = Z when NSYM = 1, and 
2MM-1 y 
= when NSYM = 2, 
MM-1 


KR(1), KR(2), and KR(3) are the axial 
stations at which the rings are located. 
(See the FORTRAN listing of Subroutine 
WIDERIL.) If there are no rings, omit 


DATA 5. 


The nondimensional magnitudes of the 
initial imperfection are placed on these 
cards. For example, if the magnitude of 
the imperfection in each mode is one 


one-thousandth of the shell thickness 
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WI (1) 


WI (2) 


WI (INMAX) 


DATA 9: 


GNU 
ALPHA 


CDAMP 


and the Fourier number n = N(L), then 
J h 


WIL(L) = TOG x a" 

The Fourier coefficient of the initial 
imperfection in the axisymmetric mode, 
n = "0, S@ Normal Wl.) 

The Fourier coefficient of the initial 
imperfection in the mode whose Fourier 
number is n = N(2). (See explanation 

of DATA 12) 

The Fourier coefficient of the initial 
imperfection in the mode whose Fourier 
number is n = N(INMAX) 

If INMAX: 


1. Is less than 15, DATA 8 must be omitted. 


2. Is less than 8, DATA 7 must be omitted. 


The nondimensional length = L/a. Sub- 
routine START sets D equal to the non- 
dimensional spacing of the finite 
difference stations. 


POLSSoOn S “Eatio 


h/a 

the nondimensional viscous damping coef- 
(1-v*) ,% 

PLeGLenc CDAMP = C i , where C is 


the viscous damping coefficient. 
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DATA 10: 


PO 


S,OMEGA 


DATA 12: 
N(1) = 0 


neCZ) 


N (INMAX) 


DATA 13: 


DT 


(As required by subroutine PRLOAD) 


- The maximum pressure 


- The decay rate 


In dimensional units, 


_ aan = sale aitees 
P(t) = Bs e = Py € 
Then 
PO =e a(l-v) 
O Eh 
ee oe 
E cer E 
blank. S and OMEGA are available to the 


user in preparing a subprogram to replace 


PRLCAD:. 


The Fourier numbers of the initial imper- 
fection. These numbers must be in one 
to one correspondence with 
WI(1),....,WI(INMAX) on DATA 6, DATA 7, 
and DATA 8. N(1) must always be equal to 


Z2E10 


The time increment to be used for the first 
time step (value used in calculations of 


Chapter V = .05) 


Zor 


DTMAX ~ The upper limit on the time increment. 
When BETA = 1/6, DTMAX should not be 
larger than one-fifth of the shortest 
period of any mode in the problem. 
Normally, the axisymmetric mode has the 
Shortest period. (value used in calcu- 


lations of Chapter V = 0.7) 


DATA 14: This data card provides the capability 
of solving several problems during a 
Single computer run. 
LAST - Ter LAST ="77,, 1f one or more™erooWen 
data sets follow. 
2. LAST = .... if the present data 


set is the last one 


Truncation of the Series 

The user prescribes which modes are in the truncated 
Fourier series by the values given to INMAX, MAXM and N(L). 
Coefficients having the same Fourier numbers as 
N (1),..-.,N(INMAX) are included in the truncated series 
for the displacements; and if MAXM is equal to INMAX, no 
other modes will be considered in the solution. If MAXM 
is greater than INMAX, subroutine MODES will add additional 
modes to the series for the displacements until the total 
number of modes iS equal to MAXM. For example, if 


INMAX = 2, MAXM = 4 and N(1l) = O, and N(2) = 3, subroutine 


Zoe 


MODES will add modes six and nine to the series for the 
displacements with zero taken as the initial imperfection 


ime meoades Six andunane. 
EV USER ALTERATIONS TO OBTAIN UNCOUPLED SOLUTIONS 


Elimination of Nonlinear Feedback to Axisymmetric Mode 

In Chapter V a procedure was outlined for identifying 
the Mathieu modes. The procedure involves computing the 
response of the flexural modes whose frequencies are close 
to one-half the frequency of the axisymmetric mode. [In 
setting up the program for such an identifying run, the 
nonlinear terms in the equations of the axisymmetric mode 
are set equal to zero. The following alterations of the 
program will eliminate the appropriate nonlinear terms. 

In subroutine LAMI1N2: 
After the statement 

oe DO? he i MNMAX® 
insert 

Gh Td 20 

100 CONTINUE 
The second inserted statement is necessary to avoid an 
error message from the Fortran compiler. 

In subroutine TNETA: 
After the statement 


9 Dd 10 L = 1,MNMAx0 


20 


insere 


GO 


TO ae 


100 CONTINUE 


The second inserted statement is necessary to avoid an 


error message from the Fortran compiler. 


Elimination of Nonlinear Coupling Between the Flexural Mdes 


In Chapter V a comparison was made between a fully 


coupled nonlinear solution and one with limited coupling. 


To eliminate the intermodal coupling between the flexural 


modes, but yet retain the coupling between each flexural 


mode and the axisymmetric mode, the following alterations 


are required: 


After the 


Insert 


After the 


Meksyeneis 


After the 


insert 


In subroutine LAMIN2: 
Staeeticmias 
IF(N(M).EQ.0) GO Td 11 
MAXL = MAXS (M) 


MAXL il 


statement 


8 MAXL = MAXD (M) 


MAXL = l 
statement 


10 IF (MAXSY(M).EQ.0) Gd Td 13 


GOMLOsIs 
100 CONTINUE 


240 


The second inserted statement is necessary to avoid an 
error message from the Fortran compiler. 
In subroutine TNETA: 
After the statements 
IF (N(M).EQ.0) Gd TO 9 
MAXL = MAXS (M) 
insert 


MAXL iL 


After the statement 
6 MAXL = MAXD (M) 
insert 
MAXL = lL 
After the statement 
8 IF (MAXSY(M).EQ.0) Gd TO 11 
insert 


COLO wr 
100 CONTINUE 


The second inserted statement is necessary to avoid an 


error message from the Fortran compiler. 
V. STORAGE REQUIREMENTS AND EXECUTION TIME 


The program requires approximately 125K bytes of 
storage space in the IBM 360 Model 67 memory when sub- 
routines AMPMAX, WI1SIN, and WI2SIN, or WIDER1 and WIDER2 
are removed from the deck. This figure includes the space 
required for the execution of the required input and out- 


put routines. 
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The approximate execution time for ycalemiating a 
completely coupled nonlinear solution until T equal two- 
hundred is shown in Figure D.2 as a function of the number 


OF modesein the Foulerer conntece 


Vi. PROGRAM Bret 


The listing of the program is contained on 


pages 244 - 279. 
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EXECUTION TIME — SECS 


4000 


3000 


2000 


l\OOO 


NUMBER OF 
FIGURE D.2 
EXECUTION TIME , ws 
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COMPUTE INITIAL RADIAL PRESSURE. 


CALL PRLCAD 
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COMPUTE INITIAL ACCELERATION. 


CALL ACELR® 
T=DT 


L 
D 
COMPUTE RADIAL PRESSURE AT END OF FIRST TIME STEP. 


CALL FRLOAD 


BEGIN FIRST ITERATION. 


1 
UPDATE RADIAL DISPLACEMENT, 


1O NIT 


AND ACCELERATION. 
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